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ABSTRACT 

We have performed a series of systematic tests to evaluate quantitatively 
the effects of spurious transport in three-dimensional smoothed particle 
hydrodynamics (SPH) calculations. Our tests investigate (i) particle diffusion, 
(ii) shock heating, (iii) numerical viscosity, and (iv) angular momentum 
transport. The effects of various program parameters on spurious mixing and 
on viscosity are investigated. The results are useful for quantifying the accuracy 
of the SPH scheme, especially for problems where shear flows or shocks are 
present, as well as for problems where true hydrodynamic mixing is relevant. 

We examine the different forms of artificial viscosity (AV) which have been 
proposed by Monaghan, by Hernquist & Katz, and by Balsara. Our tests 
suggest a single set of values for the AV parameters a and (3 (coefficients of 
the linear and quadratic terms) which are appropriate in a large number of 
situations: a ~ 0.5, f3 ~ 1 for the classical AV of Monaghan, a ~ (3 ~ 0.5 for 
the Hernquist & Katz AV, and a ~ (3 ~ 7/2 for the Balsara AV (where 7 is the 
adiabatic index). However, we also discuss how these choices should be modified 
depending on the goals of the particular application. For instance, if spurious 
particle mixing is not a concern and only weak shocks (Mach number M. ^ 2) 
are expected during a calculation, then a smaller value of a is appropriate. 
Somewhat larger values for a and (3 may be preferable if an accurate treatment 
of high Mach number shocks (Ai ^ 10) is required. We find that both the 
Hernquist & Katz and Balsara forms introduce relatively small amounts of 
numerical viscosity. Furthermore, both Monaghan's and Balsara's AV do well 
at treating shocks and at limiting the amount of spurious mixing. For these 
reasons, we endorse the Balsara AV for use in a broad range of applications. 
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1. Introduction 



Smoothed particle hydrodynamics (SPH) is a Lagrangian method introduced 
specifically to deal with astrophysical problems involving self-gravitating fluids moving 
freely in three dimensions. Pressure-gradient forces are calculated by kernel estimation, 
directly from the particle positions, rather than by finite differencing on a grid as in other 
particle methods such as PIC (the particle-in-cell method; see, e.g. |[]) or grid-based 
methods like PPM (the piecewise parabolic method; see, e.g. 0). This idea was originally 
introduced by Lucy |§ and Gingold & Monaghan ||, who applied it to the calculation 
of dynamical fission instabilities in rapidly rotating stars. Since then, a wide variety of 
astrophysical fluid dynamics problems have been tackled using SPH (see |5|, for reviews). 
In recent years, these have included planet and star formation |7], |||, ||, solar system 
formation pTOj] , supernova explosions [11], [12| , tidal disruption of stars by massive black 
holes [13 1, large-scale cosmological structure formation flU] , |l5f , galaxy formation [IB, I7| , 

The SPH method 



23, M 25 



stellar collisions [|18|, |l!| and binary coalescence [^, [H], | 
itself has also undergone major advances. Most notably, artificial viscosity (AV) has been 
incorporated [f2^, [Tj], |3"0|| , as well as powerful algorithms for the calculation of 

self-gravity including particle-mesh methods |HJ and tree algorithms |32], ^7], |33"|. 



We have performed systematic tests of the SPH method. In particular, we concentrate 
on the examination of spurious transport, including the motion of SPH particles introduced 
as a numerical artifact of the SPH scheme. Many applications require a careful tracing 
of particle positions, and in these cases it is essential that the spurious diffusion of SPH 
particles is small. For example, SPH calculations can be used to establish the amount of 
composition mixing during stellar collisions j 



18, 19], which is of primary importance 



in determining the subsequent stellar evolution of the merger remnant (see, e.g., [j3~5H ). 
Since some of the mixing observed in a SPH calculation is always spurious, the observed 
amount of mixing is an upper limit to the actual amount. Low-resolution SPH calculations 
in particular tend to be very noisy, and this noise can lead to spurious diffusion of SPH 
particles, independent of any real physical mixing of fluid elements. 

SPH particles are coupled via a smoothing kernel, and there are therefore strong 
interactions among neighboring particles. Regardless of the type of fluid being simulated, 
the physical analogue of a system of SPH particles is a crystal, liquid, or very imperfect gas 
(depending on the noise level in the calculation) but never an ideal gas of noninteracting 
particles. These particle interactions introduce a numerical viscosity into the SPH scheme 
and allow for the spurious exchange of momentum and angular momentum among shear 
layers. We have studied and measured this viscosity, both in the context of a pure shear 
flow constructed in a periodic box with slipping boundary conditions, and in a rapidly, 



-4- 



differentially rotating, self-gravitating system. 



One main goal of this paper is to present a thorough comparison of three different AV 
forms, namely those of Monaghan , Hernquist & Katz |27| , and Balsara |H| . Our tests 
of these forms include a modified version of the Riemann shock-tube problem in which 
periodic boundary conditions are imposed. For each of the AV forms, we investigate how 
the AV parameters can be adjusted to achieve an accurate description of shocks, while still 
controlling spurious mixing and shear viscosity. 

In addition, we test the effects of varying a number of SPH-specific parameters and 
schemes, including: the number of neighbors A^v, the choice of evolution equation (energy 
vs. entropy), and the type of advection algorithm. Other tests of SPH include those by 
Hernquist & Katz [27 and Steinmetz & Miiller [36|. In addition, comparisons between 



SPH and Eulerian codes have been presented in the literature in a variety of contexts: 



stellar collisions p7| , cosmology 
shock-tube tests 



, rotating stars pJ, coalescing neutron stars pE0 |, and 



OL 



Many different implementations of SPH exist (e.g. [[JT], |27|, |4"T||), and in §2 we give a 
brief description of the more popular schemes. The degree of spurious diffusion of SPH 
particles is quantified by diffusion coefficients which we measure in §3. In §4 we examine 
particle diffusion in simple self-gravitating system. In §5 we examine how well various 
strengths and forms of AV handle shocks in a modified version of the Riemann shock-tube 
test with periodic boundary conditions. In §6, we measure shear viscosity and examine the 
spurious effects of AV on the transfer of angular momentum. It is the tests of §5 and §6 
upon which we base our comparison of the various AV forms. Finally in §7 we summarize 
and discuss our major results. 



2. Numerical Method 



2.1. Density, Pressure, and Entropy 

An SPH particle can be thought of as a Lagrangian fluid element. Associated with 
particle % is its position r^, velocity v* and mass m,. In addition, each particle carries 
SPH-specific parameters including a purely numerical "smoothing length" specifying the 
local spatial resolution. An estimate of the fluid density at is calculated from the masses, 
positions, and smoothing lengths of neighboring particles as a local weighted average, 

Pi = J2m j Wi j , (1) 
j 
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where the symmetric weights = Wji can be calculated from the method of Hernquist & 
Katz ^7|, as 



W, 



1.1 



(2) 



Here W(r, h) is a smoothing (or interpolation) kernel, for which we use the second-order 
accurate form of Monaghan & Lattanzio |4lJ , 



W(r, h) 




i-fOO' + iOO 3 ; °<£<1 

2 - 



1 < I < 2, 



(3) 



f > 2. 



Depending on which evolution equation is integrated (see equations Q2DD and (^T[) 
below), particle i also carries either the physical parameter Ui, the internal energy per unit 
mass in the fluid at r$, or Aj, the entropy variable, a function of the specific entropy in the 
fluid at rj. Arbitrary equations of state (e.g. adiabatic, isothermal, even equations of state 
for metals and rocky materials; cf. p2| ) are permitted in SPH. The calculations presented 
in this paper use, unless otherwise noted, polytropic equations of state with 7 = 5/3, 
appropriate for an ideal monatomic gas. The pressure at is therefore calculated either as 



Pi = (7 - l)piUi, 



(4) 



or 



Pi = AipJ. 

We define the specific entropy of particle i to be 



1 



In 



Pi 



7 



(5) 



(6) 



and consequently the total entropy of the system S = J2i m i s i- Equation (^) is a definition 
of convenience: we refer to the quantity Sj as entropy, even though it differs from the true 
thermodynamic entropy (which depends on the composition of the fluid being represented). 
Although both Sj and the true thermodynamic entropy are conserved in adiabatic processes, 
it is Si which arises naturally when studying the dynamical stability of self-gravitating 
fluids. 



2.2. Dynamic Equations and Gravity 



- 6- 



Particle positions are updated either by 



(7) 



or the more general XSPH method 



E 



Pij 



-Wi, 



where pij = (pi + Pj)/2 and e is a constant parameter in the range < e < 1 P5| . Equation 
(|), as compared to equation (0), changes particle positions at a rate closer to the local 
smoothed velocity. The XSPH method was originally proposed as a means of decreasing 
spurious interparticle penetration across the interface of two colliding fluids. 



The velocity of particle i is updated according to 



v, = ap^ + ai SPB > 



where a\ Grav ^ is the gravitational acceleration and 



(SPH) 



Pi 



+ 



+ IT 



(9) 



(10) 



Various forms for the AV term are discussed below. The AV term ensures that correct 



jump conditions are satisfied across (smoothed) shock fronts, while the rest of equation ([10] 
represents one of many possible SPH-estimators for the acceleration due to the local 
pressure gradient (see, e.g., |43||). 



To provide reasonable accuracy, an SPH code must solve the equations of motion of a 
large number of particles (typically N >> 1000). This rules out a direct summation method 
for calculating the gravitational field of the system, unless special purpose hardware such 



as the GRAPE is used ||T7| , f44| . In most implementations of SPH, particle-mesh algorithms 
31| , |20| , |45| or tree-based algorithms |27|, [4(| are used to calculate the gravitational 
accelerations a\ Grav \ Tree-based algorithms perform better for problems involving large 
dynamic ranges in density, such as star formation and large-scale cosmological calculations. 
For problems such as stellar interactions, where density contrasts rarely exceed a factor 
~ 10 2 — 10 3 , grid-based algorithms and direct solvers are generally faster. Tree-based and 
grid-based algorithms are also used to calculate lists of nearest neighbors for each particle 



exactly as in gravitational iV-body calculations (see, e.g., J47|). 



Our SPH codes are slightly modified versions of codes originally developed by Rasio 
48| , with implementations similar to those adopted by Hernquist & Katz . Our 3D code 



has the option of including gravity, and calculates the gravitational field by a particle-mesh 
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convolution algorithm which uses a grid-based FFT solver [|7, |5 . More specifically, the 
smoothed density sets the values of the source term for Poisson's equation at grid points. 
The FFT-based convolution algorithm then solves for the gravitational potential on that 
grid. Forces at grid points are obtained by finite differencing, and then interpolated onto 
the particle positions. We have found that, for our tests involving self-gravitating fluids, it 
is relatively easy to make the gravity accurate enough that it is not a significant source of 
error. Therefore, the results of this paper can be applied to any SPH code regardless of its 
gravitational scheme. 



2.3. Artificial Viscosity 

We now present three commonly used AV forms which are tested in this paper. In §7.2 
and §7.3 we will discuss the results of these tests, while in §7.4 we discuss which of the AV 
forms performs best in which circumstances. 

For the AV, a symmetrized version of the form proposed by Monaghan is often 
adopted, 

n - = l , (li) 

Pij 

where a and j3 are constant parameters, Cjj = (cj + c 3 -)/2 (with q = (jPi/ Pi) 1 / 2 being the 
speed of sound in the fluid at r^), and 



if (v,; -v,)- (r, ; - rA < 



fj^ = 1 h^-r^/^W) ^ ^ (12) 



if (V,; - Vj) ■ (Ti - Tj) > 



with hij = (hi + hj)/2. We will refer to viscosities of this form as the "classical" AV. 
This form represents a combination of a bulk viscosity (linear in //y) and a von Neumann- 
Richtmyer viscosity (quadratic in /i^). The von Neumann- Richtmyer AV was initially 
introduced to suppress particle interpenetration in the presence of strong shocks. Morris 



& Monaghan [30] have recently implemented equation (12) with a time varying coefficient 
a, and with (3 = 2a. Our tests will demonstrate that, for constant a and (3, equation ( |TTD 
performs best when a ~ 0.5, (3 ~ 1, and rj 2 ~ 10~ 2 , although, as discussed in §7.4, these 
choices should be adjusted to fit the particular goals of an application. 

Another form for the AV, introduced by Hernquist & Katz WR calculates U.^ directly 
from the SPH estimate of the divergence of the velocity field: 

f % + % if (v« - v 7 ) • (r< - r.-) < 
^ = n Pl -f > n ' (13) 

[ if (Vj - Vj) ■ [Ti - Tj) > 
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where 

_ \apidhi\V- + f3 Pi h*\\7 ■ v|? if(V-v) i <0 
* ~ j if (V • v) 4 > li4j 

and 

(V • v), = - £ m^i - v.) ■ V<t% (15) 
Pi j 

We will refer to this form as the HK AV. Although this form provides a slightly less accurate 
description of shocks than equation QTTD, it does exhibit less shear viscosity. Our tests show 
that a ~ (3 m 0.5 is an appropriate choice for the HK AV for a broad range of circumstances 
(see §7.4). 



More recently, Balsara [2G] has proposed the AV 



n« = (^ + |) (-upij + PpI) , (16) 



where 



13 \p1 pV 

(v '~ VjHr '~ rj) &±Il if ( Vi - v,-) • ( r , - r ) ( ) 



^ = J Mlr.-r.P/^+r? 2 ) 2c ^ 1 * jV ' 1 * jV . (17) 

[ if (v< - v,) • (r< - r y ) > 

Here /j is the form function for particle i, defined by 

fi 





|V- v| 


i 


|V- v 


i + l 


V x v 


li + rfci/hi, 



where the factor rf ~ 10 4 — 10 5 prevents numerical divergences, (V • v)j is given by 
equation flT5|) , and 

(V x v)* = - Yl m i(vi ~ Vj) x ViWij. (19) 
pi j 

The function /j acts as a switch, approaching unity in regions of strong compression 
(|V • v|j >> |V x v|j) and vanishing in regions of large vorticity (|V x v|j >> |V • v|j). 
Consequently, this AV has the advantage that it is suppressed in shear layers. Throughout 
this paper we use rj' = 10~ 5 , a choice which does not significantly affect our results. 
Note that since (pi/pf +Pj/Pj) ~ 2cf,-/ (jPij), equation (|Tfi|) resembles equation flTT| ) when 
|V • v|i » |Vxv|j, provided one rescales the a and (3 in equation (ffH) to be a factor 
of 7/2 times the a and (3 in equation fllTl) . We will show that a~/5~7/2is often an 
appropriate choice for the Balsara AV. 



2.4. Thermodynamics 
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To complete the description of the fluid, either Ui or A4 is evolved according to a 
discretized version of the first law of thermodynamics: 



ditj _U (PL_L.Pi 
dt P) 



- + ^ + (vi-v^O-ViWy, (20) 



or 

<1A ' '"tE m 3 % ( v * - v j) • ViWy. (21) 
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dt 2p! 

We call equation fl2~0|) the "energy equation," while equation (|21|) is the "entropy equation." 
Which equation one should integrate depends upon the problem being treated. Each has its 
own advantages and disadvantages. For instance, thermodynamic processes such as heating 



and cooling [TJ[] and nuclear burning can be incorporated more easily into the energy 
equation. 



The derivations of equations ( P0|) and ( f21|) neglect the time variation of hi. Therefore 
if we integrate the energy equation, even in the absence of AV, the total entropy of the 
system will not be strictly conserved if the particle smoothing lengths are allowed to vary 
in time; if the entropy equation is used to evolve the system, the total entropy would then 
be strictly conserved when Ely = 0, but not the total energy [|8, |50[. For more accurate 



treatments involving time- dependent smoothing lengths, see Nelson & Papaloizou |52 
and Serna et al. 



There are many other equivalent forms of the basic SPH equations which reduce to the 
correct fluid equations in the limit N — > 00, hi — > 0. However, most of them will satisfy 
their associated conservation equations only approximately, i.e., up to errors which tend 
to zero only in this limit. In contrast, the above equations have the virtue of conserving 
energy and momentum exactly, independent of the number of particles used, as long as the 



smoothing lengths are held fixed (eg., [48]). Of course, in the numerical solution, errors will 



still be introduced by the time-integration scheme. 



2.5. Integration in Time 



For a stable time integration scheme, the timestep must satisfy a Courant-like condition 
with hi replacing the usual grid separation. For accuracy, the timestep must be a small 
enough fraction of the system's dynamical time. We calculate the timestep as 

At = C N Min(Ati, At 2 ), (22) 
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where the constant dimensionless Courant number Cm typically satisfies 0.1 < Cn ^ 0.8, 
where 

At x = Min, {K/vi) 1 ' 2 , (23) 
and where for At 2 we use one of two types of expressions, the simplest being 

= ((^W) • < 24 > 



In the presence of strong shocks, equations such as fl23j) can allow for fairly large entropy 
changes in a single timestep when Cat is large. This problem can be eliminated by using 



smaller Cat, or by adopting a more sophisticated expression introduced by Monaghan 28 



At 2 = Min^ f h ^— r > , r I • (25) 

\Ci + l.2aci + 1.2/3MaXj|//y | J v ; 



If the Hernquist & Katz AV [eq. flT3p] is used, the quantity MaXj|^-| in equation (25]) can 
be replaced by hi\V ■ v|$ if (V • v)j < 0, and by otherwise. By accounting for AV-induced 
diffusion, the a and (3 terms in the denominator of equation (p25"l) allow for a more efficient 
use of computational resources than simply using a smaller value of C^. In this paper, we 
will label the timestep routine by an S (for "simple") when we implement equations (|22|), 
(p3|), and (|24]) , and by an M (for Monaghan) when we implement (|22"D , (^3|) , and (f2~5l). 

The evolution equations are integrated using a second-order explicit leap-frog scheme. 
Such a low order scheme is appropriate because the dominate source of error for the 
evolution is the noise in particle interactions due to numerical discreteness effects. Other 
details of our implementation, as well as a number of test-bed calculations using our SPH 
code, are presented in Rasio & Shapiro |54|, p0 . 



2.6. Smoothing Lengths and Accuracy 



The size of the smoothing lengths is often chosen such that particles roughly maintain 
some predetermined number of neighbors N^. Typical values of range from about 20 to 
100. If a particle interacts with too few neighbors, then the forces on it are sporadic, a poor 
approximation to the forces on a true fluid element. In general, one finds that, for given 
physical conditions, the noise level in a calculation always decreases when is increased. 

At the other extreme, large neighbor numbers degrade the resolution by requiring 
unreasonably large smoothing lengths. However, higher accuracy is obtained in SPH 
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calculations only when both the number of particles N and the number of neighbors 
are increased, with N increasing faster than Njy so that the smoothing lengths hi decrease. 
Otherwise (e.g., if N is increased while maintaining constant) the SPH method is 
inconsistent, i.e., it converges to an unphysical limit [58]. The choice of Nn for a given 
calculation is therefore dictated by a compromise between an acceptable level of numerical 
noise and the desired spatial resolution (which is ~ h oc 1/iVjy in d dimensions) and level 
of accuracy. 



3. Simple Box Tests 



3.1. Measuring SPH Particle Diffusion 



Simulations of a homogeneous volume of gas, at rest and in the absence of gravity, 
provide a natural environment to examine spurious diffusion of SPH particles. In the ideal 
simulation of a motionless fluid, no SPH particles would move, and the thermodynamic 
variables would remain constant. However, there is always some level of noise in an SPH 
system, and this leads to the spurious motion of particles even in the absence of any bulk 
flow. 

In order to model such a system, we introduce periodic boundary conditions in a 
cubical box, adopting the standard technique of molecular dynamics (cf. [[55]]): whenever an 
SPH particle leaves the box, it is reintroduced with the same velocity vector on the opposing 
face, directly across from where it exited. Particles with smoothing kernels extending 
beyond a side of the box can have neighbors near the opposing side, once periodicity is 
taken into account. More precisely, particle j has particle i as a neighbor if there exists 
integers k, I and m such that the position (xj + kL, jji + IL, Zi + mL) is within a distance 
2hj of ( 

x ji Vji z j)i where L is the length of the box. Unless otherwise noted, the calculations 
presented in this section employ equal mass particles, all with the same time-independent 
smoothing length h chosen such that the average number of neighbors is 20, 32, 48 
or 64. The total number of particles N in the box is unimportant, as long as it is large 
enough that surface effects can be neglected. To ensure this, we always choose N such that 
L//i> 16. 

For the diffusion tests of this section, the natural units are given by n — c s — 1, 
where n is the number density of SPH particles and c s is the local sound speed. With this 
choice, velocities are in units of c s , distances are in units of n -1 / 3 , and times are in units of 
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n~ l ^c~ l . In practice, we implement c s = 1 by choosing the entropy variable A = p 1_7 /7. 
Furthermore, the mass of the particles is chosen such that the cubical box contains unit 
mass: M = Nm = 1. Since the local number density and sound speed are known in any 
SPH calculation, these units make our results applicable to many contexts. 

After positioning the particles on a regular lattice and assigning their velocities (with 
zero net momentum), we allow the system to evolve, without AV. Although each SPH 
particle represents a physical fluid element with a certain temperature and density, the 
SPH particles themselves have their own numerical "temperature" (due to the particle 
velocity dispersion) and number density. While there is an obvious correlation between the 
number density of the SPH particles and the density of the gas being represented, no such 
correlation exists between the numerical temperature of the SPH particles and the physical 
temperature of the gas being simulated. Regardless of the initial velocity distribution 
chosen, the velocities ultimately settle into an equilibrium Maxwell-Boltzmann distribution 
(see Figure [l]) with some root mean square particle velocity v rms which quantifies the 
noise level, or numerical temperature, of the system. The energy exchange which causes 
this thermalization is due to the strong coupling between neighboring particles through 
equation ([U]). We have also found that the velocity distribution in real calculations tends 
to be roughly a Maxwellian centered on the local smoothed velocity. 

The level of diffusion is quantified as follows. Once the velocity distribution has settled 
into an equilibrium Maxwellian, we record the positions of all particles. Since ideally the 
particles would not move far from their initial positions, it is then easy to monitor the 
mean square spurious diffusion distance 5 2 as a function of time t (properly accounting 
for particles which cross the faces of the box). At late times the mean square deviation 
5 2 increases at an nearly constant rate, so that the system obeys the usual diffusion 
equation 5 2 = Dt, and the diffusion coefficient D = d5 2 /dt, evaluated at late times, is easily 
measured. (In molecular dynamics, the diffusion coefficient D is sometimes defined to be a 
factor of six smaller than in our definition.) As an example, Figure |2| shows 5 2 and d5 2 /dt 
for a system with an equilibrium v rms = 0.069; it is clear that dd 2 /dt is essentially constant 
at late times, and we measure D « 0.024. 

Figure |3] shows the diffusion coefficients D for various v rms and for Nn =20, 32, 48 and 
64. Not surprisingly, spurious diffusion increases as v rms increases. Note that, for a given 
N^i there is a critical noise level below which the diffusion coefficient D is essentially zero. 
In this regime, the SPH particles settle into a regular lattice and oscillate around their 
equilibrium positions, and we say the system has "crystallized" (see §3.2). There seems 
to be a crystallization point for all the curves at some critical velocity dispersion v cr > 0. 
The trend is for v cr to decrease as Nn increases. During the dynamical phase of real 
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applications, AV typically keeps the noise level low enough that the numerical temperature 
is at most slightly above that required for crystallization. 

The diffusion coefficient is not always a unique function of Nn and v rms , but can 
also depend on the history of the SPH particles. To demonstrate this we started the 
particles on various types of lattices. Figure |] shows the measured values of the diffusion 
coefficient D in the crystallization regime for systems of particles which began in either face 
centered cubic (dashed lines) or a simple cubic (solid lines) configurations. There is a clear 
dependence on the system's history in this regime, making it impossible to define a precise 
crystallization velocity dispersion. Note that all of the data points in Figure |] have a small 
diffusion coefficient, D < 0.025. Well above the crystallization noise level (that is, outside 
of the region displayed in Figure 4) the diffusion coefficient is largely independent of initial 
conditions; that is, there is negligible history dependence for sufficiently large v rms . 

The diffusion coefficients shown in Figures |3] and [| are measured while integrating the 
entropy equation (pl|) with a Courant number = 0.4 and with the S timestep algorithm 
[see eqs. (0), (0), and (0)]. However, measurements which use the energy equation (|20| ) 
or different Courant numbers, or both, give similar coefficients, provided only that the 
Courant number is small enough that the integration routine is stable. 



3.2. Lattices of SPH particles 



By experimenting with various lattice types as initial conditions in the simple box 
tests, we have found that not all equilibrium configurations of SPH particles are stable. 
For example, simple cubic lattice configurations are unstable to perturbations, while other 
lattice types, such as hexagonal close-packed, are stable. If the particles begin motionless 
and slightly perturbed from equilibrium simple cubic lattice sites, they achieve a non-zero 
noise level and readjust their positions to a different, preferred lattice type (see Figure [5]). 
Although the instability develops more slowly for smaller CV, it cannot be avoided 
altogether. 

For a few of our simple box tests, we allowed the smoothing lengths hi to vary both in 
time and in space, without including the corrections in the evolution equations described 



by Nelson & Papaloizou [51, 52] and Serna et al. [53 1 . The system's behavior is greatly 
affected: there is a secular, spurious increase in the total energy E. Almost all of this 
spurious energy is kinetic. If the AV is active during such runs, energy conservation is much 
better; however, the error then emerges as a spurious entropy increase (see Figure |6|). The 
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AV run in Figure y used a = 1, (3 = 2, rf = 0.01 and the classical form of AV; both runs use 
Cn = 0.8 and an initial Maxwell-Boltzmann velocity distribution with a velocity dispersion 
v rms = 0.107. 



In many SPH applications, shocks play an important role in the dynamics. Therefore, 
understanding how various AV schemes affect the level of spurious diffusion is essential. 
A uniform SPH gas is not an appropriate arena to study this effect, since the AV quickly 
solidifies the particles into a lattice structure. In a calculation with AV but without shocks 
or shear, the diffusion coefficient D is always essentially zero (see Figures [j] and ^j, since 
diffusion occurs only transient. 

We can derive approximate analytic expressions for the artificial viscous dissipation 
timescale by dimensional analysis on the AV term in equation fllQ|) . Here we focus on the 
classical AV [eq. pf ; in §6.2 we will analyze all three AV forms in a different context. 
Beginning with equation (|T2"D, we note that since |rj — rj| ~ hij we have fiij ~ Av, where 
Av is a typical relative velocity of neighboring particles. If, in the vicinity of particles i and 
j, the sound speed is c s and the density is p, then equation ( |Tl"D gives us ~ — a Avc s /p 
if f3Av << ac s (as is typically the case in the absence of shocks). If the local number 
density of particles is n, then a typical particle mass mj ~ p/n, and iVjW 7 ^'! ~ n/(hNN). 
Combining these expressions, we find that the acceleration of particle i due to the AV is 



■AV 



ac* Av 



v 



= I~E m^ViWijl ~ TT^f, (26) 



hN N 



where we have assumed that the sum over Nn terms in equation (|Tt]) scales as N^ 2 since 
there is no preferred direction for ViWij. 

The artificial viscous dissipation timescale r is then just v/v AV , where v is a typical 
particle velocity. For the simple box tests we have v ~ Av ~ v rms , so that the viscous 
timescale is 

h N l/2 / ^ \ V3 /v 5 / 6 
T „ = f A\ £h^ n -^ c -\ (27) 



ac s V327T/ a 

Our numerical results agree well with this simple expression. For a = 1 and Nn = 32, 
equation (B71) gives a timescale r ~ 6n _1 / 3 c~ 1 , which is in reasonable agreement with the 



time it takes to form a lattice (i.e. the timescale on which the kinetic energy drops to zero) 
in the case presented in Figure 0. Although the timescale depends on both and the AV, 
it is always quite short: typically just a few sound crossing times between neighboring SPH 
particles. 
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4. Polytrope Tests 



Applications of SPH often involve self-gravitating systems with significant density 
gradients. The results of our simple box tests can be applied to such calculations, which 
we will demonstrate by considering a set of equilibrium n = 1.5 polytropes (spherical 
hydrostatic equilibrium configurations with p = const x p 1+1 / n ) all with mass M and radius 
R, but modeled with various total numbers of particles N and neighbor numbers N^. In 
this section, all calculations implement the simple timestep routine given by equations 
([22l)-(|24l) and have no AV. The natural units are given by G = M = R = 1, so that 
consequently the unit of time is (R 3 /GMY^ 2 . 

We relax the polytrope to equilibrium by applying an artificial drag force which opposes 
motion for 20 time units. We then remove the drag force and record the particle positions. 
Ideally, the particles would remain stationary. However, as expected from the results of 
§3.1, these particles spuriously diffuse from their starting positions, and this diffusion is 
easy to monitor. By periodically noting the particle velocity dispersion v rms , we can apply 
the simple box test results to get an "instantaneous" value for the diffusion coefficient D by 
interpolating between data points in Figure |3[ In this way, we "predict" the mean square 
displacement S 2 by a simple, numerically evaluated integral, 



and then compare this prediction to the actual, measured mean square displacement. 

Figure ^ shows, as a function of time, the mean square spurious displacement for the 
innermost 6400 particles in an n = 1.5 polytrope modeled with N = 13949 equal mass 
particles, each with N N = 48 neighbors on average. We do not track the particles of the 
outer layers here, since they are subject to an effect which we do not attempt to model: 
when such a particle diffuses outward beyond the surface, gravity pulls it back, making the 
actual diffusion distance somewhat smaller than predicted. For those particles which always 
remain inside the surface, gravity is everywhere balanced by pressure gradient forces, so 
that the rate of diffusion is essentially the same as in our simple box tests. The usual 
advection scheme equation (|7]) was used for the calculation presented in the top frame of 
Figure |9|, while the XSPH equation (||) with e = 0.5 was used in the bottom frame. The 
"predicted" mean square displacement (dashed curve), as calculated from equation (^8|), 
agrees well with the actual square displacement (solid curve). To obtain the predicted curve 
in the XSPH calculation, the root mean square of the right hand side of equation (|8|) was 
used in place of v rms when determining the diffusion coefficient D. The Courant number 
CV = 0.8 and the simple timestep routine determine the integration timesteps for both 
cases. 




(28) 
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The slight differences between predicted and actual displacements arise because of 
our interpolating to obtain D and because our diffusion coefficients are only approximate 
in the crystallization regime (due to history dependence). Since the SPH particles are 
melting out of their crystalline phase around t w 10, our values for D are overestimated 
then. The XSPH advection method does indeed diminish the amount of spurious diffusion: 
the final (t = 30) mean square displacement for the XSPH calculation is nearly one fourth 
of the value from the simple advection scheme. However, one must be careful when using 
XSPH: using too large of an e can cause certain modes to become numerically unstable. 
For instance, for the extreme case of e = 1 we are not able to evolve an equilibrium n = 1.5 
polytrope without the integration becoming unstable. 

In order to test the importance of the Courant number Cn, we evolved a set of 
equilibrium n = 1.5 polytropes using several values of Cn between 0.1 and 1.6. In all cases 
we turned off the AV, used iV = 13949 equal mass particles each with Nn = 48 neighbors 
on average, allowed the smoothing lengths to vary in space and time, used the "simple" 
timestep routine, and monitored three measures of error: the fractional (spurious) change 
in total energy AE/E, the velocity dispersion (v/c s ) rms , and the mean square diffusion 
distance 5 2 . As Cn increased from 0.1 to 1.1, these errors, evaluated at t — 25, increased 
only very slightly: AE/E increased from 0.014 to 0.017, (v/c s ) rms from 0.13 to 0.14, and 
5 2 /R 2 from 0.02 to 0.03. For Cn = 1-2, the integration becomes unstable, with the errors 
at t = 25 then being AE/E = 0.7, {v/c s ) rms = 0.3, and 5 2 /R 2 = 0.2. This result suggests 
that in certain cases, for fixed computational resources, it may be more efficient to use 
a relatively large Courant number like Cn = 0.8 and more particles, rather than a small 
Courant number like Cn = 0.3 and fewer particles. 

Figure |H] shows AE/E, (v/c s ) rms , and S 2 /R 2 at t = 25 for a set of calculations with 
Cn = 0.8 and various N N - Here the n = 1.5 polytropes are modeled by either N = 30000 
particles (circular data points) or N = 13949 particles (square data points). For a given 
Nn, the N = 30000 models always have larger accumulated errors: as N is increased, 
one must also increase Nn in order for the SPH scheme to remain accurate. Although 
increasingly larger Nn results in increasingly smaller errors, this does not mean one should 
strive to use as large a value for Nn as possible. Large Nn yields large smoothing lengths 
and hence poor spatial resolution. The optimal Nn must be determined by a compromise 
between the competing factors of accuracy and resolution, and depends on the particular 
application. Nevertheless, we can place very loose constraints on how fast the optimal Nn 
should be increased as N is increased. From Figure |TU| we see that in going from iV = 13949 
to iV = 30000 we need to increase Nn by at least (very roughly) 15% in order to prevent 
the errors from increasing. This corresponds to a scaling Nn oc N q with 0.2 <^ q < 1, 
assuming a power-law relation. The upper limit of 1 on q stems from the requirement that 
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the smoothing lengths must decrease as N and increase. 

We have performed our diffusion tests using equal mass particles. Sometimes, however, 
SPH calculations use particles of unequal mass so that less dense regions can still be 
highly resolved. Unfortunately, the more massive particles tend to diffuse to the bottom of 
the gravitational potential more so than less massive ones. In other words, each particle 
has a preferred direction to diffuse, and in a dynamical application this direction can be 
continually changing. As an example, we evolved an equilibrium n = 1.5 polytrope in which 
the SPH particles initially in the envelope were, on average, heavier than those in the core. 
Over the course of the calculation, the heavier particles settled to the core while the lighter 
particles tended to the envelope (see Figure |TT|). Such behavior makes spurious diffusion 
more difficult to predict for calculations which use unequal mass particles. 



5. Periodic Shock- Tube Tests 



Since the simple box tests of §3 are helpful only for calculations without AV, we turn 
now to a periodic version of the ID Riemann shock-tube problem of Sod [56], a standard 



test of hydrodynamic codes and AV schemes containing many of the same qualitative 
features as real applications which involve shocks. The physical setup is as follows. 

Initially, fluid slabs with constant (and alternating) density p and pressure p are 
separated by an infinite number of planar, parallel, and equally spaced interfaces. If we 
define the unit of length to be twice the distance between adjacent interfaces, and if we let 
the x = plane coincide with one of these interfaces, then 

P = Pi, P = Pi if -| < x < 

p — p r , p — p r if < X < ~, 

where pi, pi, p r and p r are constants specifying the density and pressure of the slabs 
to the "left" and "right" of x = 0. Pressures and densities for \x\ > \ are given by 
repeatedly stacking the thermodynamic slabs described by equation (2^) along the x-axis 
to infinity, hence the name periodic shock-tube tests. At t — the interfaces are removed 
and, if pi ^ p r , a shock wave moves from the high pressure material into the low. A 
rarefaction wave also originates at each interface, propagating in the direction opposite to 
its corresponding shock. Before the initial collision of shock waves from adjacent interfaces, 
regions of five different thermodynamic states coexist and the entropy of the fluid increases 
linearly with time. A quasi-analytic solution can be constructed for these early times using 



standard methods (see, e.g., [[5^]]) and is presented in detail by Rasio & Shapiro |54 . 
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5.1. Low Mach Number Cases 



For the first set of shock-tube calculations we consider, the fluid slab to the left 
of the interface at x = initially has density pi = 1.0 and pressure pi = 1.0, while on 
the right p r = 0.25 and p r = 5/2 16 / 3 = 0.12402. Consequently this box contains 0.625 
units of mass: 0.5 on the left and 0.125 on the right. An adiabatic equation of state 
is used with 7 = 5/3, so that the entropy variable A equals 1.0 on the left and 1.25 
on the right. From equation (|6|), the initial entropy of each of the periodic cells is thus 
S = 1.5[0.51n(1.5) + 0.125 ln(1.5 x 1.25)] = 0.4220. For these initial conditions, the initial 
shock waves have a relatively low Mach number Ai »s 1.6. In these units, the speed of 
sound in the initial left hand slab is c l s = {^Pi/ pi) 1 ^ 2 = 7 1//2 , and the unit of time is therefore 
7 1//2 L/c^, where L is the length of a periodic cell (our unit of length). 

Employing AV with the form of equation (|TT]) , we obtained a good representation of 
the shock with our ID code by using the AV parameters a = (3 = 1 and rj 2 = 0.05 in the 
classical AV. The smoothing length h of the N = 2500 equal mass particles was constant 
and chosen such that the particles would have Nn = 16 neighbors on average. Our ID code 
integrates the energy equation, and uses the Monaghan timestep routine with Cn = 0.2. 
Figure [TJ| shows the density and velocity profiles as given by the quasi-analytic solution 
(solid curve) and our 1-dimensional code (dotted curve) at a time t = 0.15. As expected, 
discontinuities are smoothed over a few smoothing lengths. Figure [T3] shows the entropy [see 
eq. (P)] given by our ID SPH code (dotted curve), which nearly matches the quasi-analytic 
solution (solid curve). 

The above calculation helps establish the accuracy of our ID code, but does not 
realistically assess the accuracy of a 3D calculation, where the much smaller number of 
particles per dimension leads to a reduced spatial resolution. Furthermore, many sources 
of numerical errors, including spurious mixing, are artificially reduced when motion with 
only one degree of freedom is allowed. We test our 3D code by using it to simulate exactly 
the same physical problem: at t — 0, slabs of fluid with alternating thermodynamic states 
are separated by equally spaced planar interfaces perpendicular to the x-axis. Periodic 
boundary conditions are imposed on all six sides of a cube with faces atx = ±|,y = ±|, 
and z — ±5. We considered cases only with a constant smoothing length h « 1, and, 
unless otherwise stated, we integrate the entropy equation. 

Our calculations with the 3D code use N = 10 4 equal-mass particles. All the particles 
initially in the left hand slab have the same smoothing length, smaller than the smoothing 
length common to particles initially in the right hand slab. These smoothing lengths are 
not allowed to vary with time, and are chosen such that particles which are farther than 2h 
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from an interface have = 64 neighbors on average. Within each constant density slab, 
the SPH particles start in a stable lattice with a randomly chosen orientation (choosing the 
lattice face to be parallel to the interface would be too artificial of a setup). The initial 
conditions for each slab are constructed by randomly distributing particles in a periodic box 
of dimensions | x 1 x 1 and then slowly relaxing the system with an artificial drag force. The 
resulting lattices are preferred to initially randomly distributed particles, since a random 
distribution would introduce a high noise level not representative of real applications. 

We determine the accuracy of our calculations with the 3D code by comparing its 
results against those of the much more accurate ID code. Such 3D calculations are a 
useful and realistic way to calibrate spurious transport in simulations with AV, since the 
test problem, which includes shocks and some large fluid motions, has many of the same 
properties as real astrophysical problems. In fact, the recoil shocks in stellar collisions do 
tend to be nearly planar, so that even the ID geometry of the shock fronts is realistic. The 
periodic boundary conditions play the role of gravity in the sense that they prevent the gas 
from expanding to infinity. 

Figure |T^ shows the pressure P, entropy variable A, density p and velocity v x as given 
by our ID code (solid curve) and by our 3D code (dots) at the relatively late time t = 1. 
Here the 3D calculation implements the classical AV with a = 0.5 and (3 = 1. The bar 
in the lower left corner of the uppermost frame displays the average region of influence 
(i.e. the mean diameter of the smoothing kernels) for the particles in the 3D calculation: 
the total length of this bar is 4(/i), where (h) = 0.058 is the average smoothing length. The 
3D calculation does well at reproducing the major features in the thermodynamic profiles, 
but, not surprisingly, smoothes out any small scale structure which occur on lengths scales 
shorter than a few smoothing lengths. In the regions near x = 0.1 and x = 0.4, where the 
fluid is being shock-heated, the pressure, entropy variable and density in the 3D calculation 
are double-valued due to the shock front not remaining perfectly planar throughout the 
calculation. 

Since the fluid motion in these calculations should be solely in the x-direction, spurious 
motion in the y- and ^-directions is easy to measure. Spurious motion in the x-direction 
can be studied by the following method, based on the idea that planes of fluid should 
not cross in one-dimension. That is, the shape of a composition profile should remain 
unchanged throughout a calculation. Once the shock-tube system has reached a steady 
state, we examine the distribution of the Lagrangian labels Xi(t = 0) as a function of m(x), 
the amount of mass between the interface (contact discontinuity) and x. Deviations from 
the initial profile must be spurious, so we can immediately calculate spurious displacements 
in the x-direction for individual particles. Diffusion measurements in each of the three 
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directions give similar results. 

We have studied the quality of the 3D code's results for various AV parameters and 
forms. We have completed a number of shock-tube tests which began with the same initial 
conditions described above, but with different values of the AV parameters. Varying rf 
by a factor of 25 between 0.002 and 0.050 makes little difference in the results, and we 
therefore concentrate on the effects of a and (3. All calculations described in this section 
have rj 2 = 0.01. 



Figure ITS] shows the dependence of the solution on a and (3 for the classical AV by 



plotting, as a function of time, the mean square spurious displacement in the directions 
perpendicular to the bulk fluid motion (in units of n" 1 / 3 , where n is the SPH particle 
number density), the internal energy U and the entropy S. The solid line results from our 



accurate calculation of the shock-tube problem with the ID code. In frame (a) of Figure [15 
a = while j3 is varied. In (b) (3 — and a is varied. Finally in (c) (3 — 1 and a is varied. 
Runs with a = or (3 = are interesting since they represent an AV which is either purely 
quadratic (von Neumann- Richtmyer viscosity) or linear (bulk viscosity) in /j,^, respectively, 
and these two types of AV generate different numerical viscosities (see §6). 

Table I summarizes all of our low Mach number 3D shock-tube calculations and reports 
how well each does matching the ID solution. All the calculations in Table I employed 10 4 
particles and a fixed smoothing length chosen such that the number of neighbors = 64 
on average. In Column 1, we identify the type of AV used: C for the classical AV [eq. (JTTJ)] , 
HK for the Hernquist & Katz AV [eq. (|i"3D1, and B for the Balsara AV [eq. ([16])]. Columns 
2 and list the AV parameters a and (3 (unless otherwise noted rj 2 = 0.01). Column 4 gives 
the type of timestep routine used: S for simple [eq. J23p] and M for Monaghan [eq. (EBI)] . 
Column 5 gives the Courant number C^. Columns 6 and 7 give the number of iterations 
required to reach t = 1 and t = 4, respectively. Column 8 gives the fractional deviation 
in the total energy away from its initial value: AE/E = \E(t = 4) — E(t = 0)\/E(t = 0). 
The t = 4 value of 5 2 + 5 2 , the spurious displacement squared in the direction perpendicular 
to the bulk fluid flow, averaged over all particles, is listed in Column 9. Columns 10 and 
11 give the maximum deviation in U/E and S, respectively, from that of the ID code: 
A(U/E) max = M&x\U 3D /E 3D - Ui D /E 1D \ and AS max = Max\S 3D - Si D \. 

Figure [T^ shows, as a function of time, the average square displacement perpendicular 
to the bulk fluid flow 5 2 + 5 2 , the ratio of internal to total energy U/E, and the entropy S 
for three calculations with different forms of AV: the classical AV with a = 0.5, /3 — 1 (long 
dashed curve), the HK AV with a = (3 = 0.5 (short dashed curve), and the Balsara AV with 
a — (3 = 7/2 (dotted curve). In the bottom two frames, the solid curve corresponds to our 
ID SPH code. As we will discuss in §7.4, these choices for a and (3 are our recommended 
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values. We see that all three AV forms can handle the shocks with roughly the same degree 
of accuracy, although the HK AV does allow slightly more spurious mixing and does not 
match the ID code's U/E curve quite as well. 

We also ran several low Mach number calculations with the energy equation being 
integrated. Table II compares these runs against the corresponding calculations in which 
the entropy equation was integrated. For given values of a, (3 and 7] 2 , the two schemes do 
equally well at conserving energy, at controlling particle diffusion, and at matching the 
time evolution of U/E from the ID calculation. However, integrating the energy equation 
does allow slightly larger errors in the evolution of entropy, with AS max being 0.005 to 
0.007 larger than when the entropy equation is integrated. This larger error in the entropy 
accumulates mostly at early times when the shocks are strongest. 



5.2. High Mach Number Cases 

Since many astrophysical situations involve shocks which are stronger than the low 
Mach number situation described in the previous section, we repeated shock-tube tests with 
a larger difference in pressure between the alternating fluid slabs. In particular, we initially 
set pi = 1.0 , pi = 1.0 and p r = 0.25. but reduced the pressure of the right-hand fluid slab 
to p r = 1.2402 x 10~ 3 , a factor of 100 less than in the low Mach number cases of §5.1. This 
increases the Mach number of the initial shock waves to M. ~ 13.2. The initial entropy of 
each of the periodic cells is S = 1.5[0.5 ln(1.5)] + 0.125 ln(1.5 x 0.0125)] = -0.4415. 

For our ID code, we continued to use the classical AV [see eq. (|H])] with parameters 
a = (3 = 1 and rf = 0.05. We used 2500 particles and constant (in time) smoothing lengths 
hi, chosen such that the particles have 16 neighbors initially. Figure [17] shows a comparison 
between our ID SPH code (dotted curve) and the quasi-analytic solution (solid curve) at 
a time t = 0.15. As expected, the ID code does smooth out discontinuities in the density 
over a width of a few smoothing lengths. However, the agreement between the ID code and 
the quasi- analytic solution is still very good. 

As in the low Mach number case, we can compare the results from the 3D code to that 
of the ID code, in order to evaluate the amount of spurious mixing and to determine the 
acceptable range of values for the AV parameters for our 3D calculations. Table III is the 
high Mach number equivalent of Table I. These 3D calculations employ iV = 10 4 particles 
each with Nn = 64 neighbors, as in the 3D low Mach number calculations. 

In Figure |1^ we present the results of 3D high Mach number shock-tube calculations 
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for various a and (3 with the classical AV. For all the 3D calculations in this figure, we 
chose rf = 0.01 and used the Monaghan timestep routine with Cn = 0.8. The solid line is 
the result of the ID calculation. It is apparent that the spurious displacement is smaller 
for stronger AV, as expected and as in the low Mach number tests. As also seen in the 
low Mach number tests, the case with the lowest spurious mixing (a = 5, (3 = 0) has the 
worst fit to the energy curve of the ID calculation. The entropy curve from the ID case lies 
between the cases with the high values of a or (3, and those with the low values. Therefore, 
the best choice of AV parameters will depend on the particular situation which is to be 
modeled. If spurious mixing is important to control, then a strong viscosity is favorable. 
On the other hand, if spurious mixing is not an issue, one could use a weaker AV to more 
accurately determine the evolution of the system. 

Figure [H5| shows, as a function of time, the average square displacement perpendicular 
to the bulk fluid flow 5y + 5%, the ratio of internal to total energy U/E, and the entropy S 
for three calculations with different forms of AV: the classical AV with a = 0.5, /3 — 1 (long 
dashed curve), the HK AV with a = (3 = 0.5 (short dashed curve), and the Balsara AV with 
a = (3 = 7/2 (dotted curve). In the bottom two frames, the solid curve corresponds to our 
ID SPH code. As will be discussed in §7.4, these choices for a and (3 are our recommended 
values. We see that the HK AV does allow slightly more spurious mixing and does not quite 
match the ID code's U/E curve as well. Nevertheless, all three AV forms can adequately 
treat the strong shocks of this system. 



5.3. High Mach Number Cases with 7 = 3 

Of course, not all fluids are well-described by the ideal gas (7 = 5/3) approximation. 
For example, neutron star matter is best represented by a stiff equation of state with 
7 Pd 2-3, while an isothermal gas can be described with 7 = 1. Changing the value of 
7 changes the thermodynamic properties of the material we model with SPH, which in 
turn affects the way the AV behaves. Therefore, to investigate the dependence on 7 of the 
'optimal' AV parameters, we have performed some shock-tube calculations with 7 = 3. The 
fluid slabs were set up to have the same Mach number as the previous high Mach number 
ideal gas tests (M = 13.2): p t — l,pi — l,p r = 0.25, and p r = 8.78 x 10~ 7 . The initial 
entropy of each periodic cell is S = 0.5 [0.5 ln(0.5) + 0.125 ln(0.0562/2)] = -0.3965. 

For the corresponding calculation with the 1-D code, we used the classical AV scheme 
with a — (3 — 1 and r/ 1 = 0.05, 2500 particles and 16 initial neighbors, as in the previous 
high Mach number case. For our 3-D calculations, we used 10 4 particles with 64 initial 
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neighbors, and a variety of AV parameters with all three AV schemes. We used the 
Monaghan timestep routine with Cjy = 0.3. The different 3-D calculations and their 
comparison to the ID runs are given in Table IV, and a selection of the results are shown 
in Figure 

As in the ideal gas case, spurious diffusion is smaller for stronger artificial viscosities. 
The calculations with small a show additional "wiggles" in the energy curve (see Fig. 
pOD and larger errors in energy conservation (see Table IV), suggesting the appearance of 
numerical instabilities for strong shocks treated by weak AV forms. In general, we find that 
the level of energy conservation is worse in our 7 = 3 calculations than in our 7 = 5/3 
calculations (compare Tables 3 and 4). 



6. Shear Flows 

6.1. Periodic Box Tests 

In order to model a shear flow of infinite extent, we return to a cubical box with a side 
length L = 1 and periodic boundary conditions. The boundary conditions on the x = ±|f 
and z = ±~ faces are identical to the periodic boundary conditions in the simple box tests 
of §3: when a particle crosses one of these faces it is reinserted with the same velocity at the 
corresponding position on the opposing face. On the y = ±|f faces, however, we implement 
"slipping" boundary conditions in order to maintain a velocity field with a shear flow: if 
a particle crosses a face with a velocity (v x ,v y ,v z ) at a position (x, ±~,z), it is reinserted 
with a new velocity (v x =F ^o, v y , v z ) at the position (x =F Vot + kL, =f-|, z), where t is the time 
elapsed since the beginning of the calculation and k is the integer which places the particle 
in the central periodic cell. The resulting "stationary Couette flow" has a velocity field 
close to (v y/L,0,0) (see Fig. gl]). 

Neighbor searching across the x = ±-| and z = ±|- faces is done exactly as in §3.1. 
Across the y = ±~ faces, the slipping boundary conditions are taken into account: the 
criterion for particle j having particle i as a neighbor is that there exists integers k, I and m 
such that the position (xi + kL + Ivot, yi + IL, Zi + mL) is within a distance 2hj of (xj, yj, Zj). 
In addition, the relative velocity of particles interacting across the y — ±4 boundaries is 
adjusted by vq when computing the AV term n^. In this way, particle interactions across 
the boundaries behave identically to interactions within the box. 
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Our units of distance and mass are the length of the box and the total mass in the box: 
L = 1 and M = Nm = 1, where N is the number of particles. We set the entropy variable 
A = 1 for all the particles initially. Consequently the unit of velocity in our calculations is 
7~ 1 / 2 c s , where c s is the initial sound speed, and the unit of time is 7 1//2 L/c s . 



Figure |22| shows the spurious square displacement, energies, and entropy as a function 
of time in three calculations with N = 1000, = 64, vq = 0.l7 -1//2 c s , and various forms 
of AV. The system was relaxed for the first 10 time units (without AV) towards a situation 
with (v x ,v y ,v z ) = (voy/L, 0, 0), while from t = 10 to 50 the system evolves freely with the 
slipping boundary conditions and AV. 

Notice the increase in energy once the relaxational damping is turned off: roughly a 
1% increase in E per time unit. This increase results from the slipping boundary conditions 
and, for a given AV form and AV parameters, is nearly independent of the resolution. 
Since we are moving the boundary surfaces by hand and since there is viscosity, there is 
a shear stress at the boundaries and work is being done on the system. This behavior is 
analogous to that of a truly viscous fluid forced to undergo shear flow between close moving 
boundaries (as in a viscosimeter) : the added energy goes into heating the fluid. 

Since the faces of our cubical box have surface area L 2 , the viscous force F x acting on 
the fluid inside of the y = ±L/2 faces is 

F x = ±V^L 2 = ± V v L, (30) 
dy 

where rj is the dynamic viscosity (not to be confused with the AV parameter rj 2 ). The rate 
of energy change of the system is therefore 

dE 

-fa = [ F xVx} y= _ L/2 + [F x v x } y=+L/2 (31) 
= VVoL, (32) 

Measuring the rate of energy increase therefore allows us to numerically determine the 
dynamic viscosity. This procedure for measuring viscosity is also implemented in molecular 
dynamics (eg. ||58|| ). To calculate the kinematic viscosity v from the dynamic viscosity rj, 
one simple uses u = i]/p = r]N/(Mn), where n is the number density of particles. 

In the absence of any spurious motion, SPH particles should maintain the same 
spatial coordinates y and z throughout the calculation. By monitoring motion in these 
two dimensions, we can therefore easily quantify the extent of spurious diffusion. As in §3, 
the square displacement increases linearly with time at late times. Here we measure the 
diffusion coefficient D by fitting the relation 3(5 2 + 5 2 )/2 = Dt. In practice, we determine rj 
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and D from the average slope of the energy and square displacement curves, respectively, 
between times t = 12 and t = 50. Tables V and VI list the results of a set of runs at two 
different shear velocities with = 64. We vary the AV scheme and the AV parameters, 
and monitor the time averaged velocity dispersion ((v 2 + v 2 z )/c 2 s ) between t — 12 and 50. 
We also list the viscosity r\ [as determined from eq. (H)], the diffusion coefficient D, and 
the product r/D for each case (all converted into units M = c s = n = 1 to keep our results 
applicable to general situations). In the last three columns, the number in parentheses 
"()" is the error in the last digit, or last two digits, that is quoted. The uncertainties for 
the viscosity rj and the diffusion coefficient D are determined from the root mean square 
deviation of E(t) and of Sy(t) + S 2 (t) from the best-fit linear curve. In Tables VII we 
present results from a handful of calculations with various neighbor numbers N^. All of 
the calculations use constant smoothing lengths, as well as a constant timestep dt = 0.01 so 
that fixed computational resources are available. 



6.2. Rapidly rotating, self-gravitating fluids 

Rotation plays an important role in many hydrodynamic processes in astrophysics. For 
instance, the collision of two stars typically results in a rapidly and differentially rotating 
merger remnant. Even in the absence of shocks, AV tends to damp away differential 
rotation due to the relative velocity of neighboring particles at slightly different radii. Many 
systems are best modeled as a perfect fluid, ideally with a viscous timescale r = oo. In 
such cases, any viscosity introduced by the SPH scheme is spurious. In this section, we 
consider a differentially rotating, self-gravitating fluid and analytically estimate the viscous 
timescale for each of the three AV forms examined in this paper. Our analytic estimates 
are then compared against numerical determinations of the viscous timescale. The larger 
the viscous timescale, the more closely the calculation is treating the gas as a perfect fluid. 

As our concrete example, we consider an axisymmetric equilibrium configuration 
rotating with an angular velocity Q oc xzj~ x , where the cylindrical radius w is the distance 
to the rotation axis. In this case, the magnitude of the quantity (v, — v 3 -) • (r, — r^) which 
appears in equation (|P2| ) is ~ h Av for two neighboring particles separated by ~ h, a typical 
smoothing length, where the shear velocity Av = XQh. Note that Av = for the special 
case of rigid rotation (A = 0). 

If the AV is of the form of equation ([LI]) with (3 = 0, equation ( p6|) gives v AV and 
the viscous dissipation timescale r = v/v AV = Vtw/v AV ~ wNlj 2 j '(a\c s ). Note that this 
timescale r is not directly dependent on N: increasing N while keeping fixed would 
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not therefore affect the viscous timescale in this case. For general a and (3, we analytically 
estimate from equation ( |TT] ) that 



n*i « -ji^^ - j 2 ^^, (Classical AV) (33) 
P P 

where ji and ji are dimensionless coefficients of order unity. In this case equation (f26|) 
must be replaced by v AV ~ k\ac s Av / (hN^J 2 ) + k2/3Av 2 /(hNlf 2 ), and the viscous timescale 
r = v/v AV is then given by 



_ ■ , , aAvc s [3Av 2 \ 1 / a\c s flXAv ' . 



(34) 

where fci and fc 2 are dimensionless coefficients of order unity. The ratio of the two terms 
on the right hand side of equation ( [51D tells us that the von Neumann-Richtmyer viscosity 
(corresponding to the term with f3) yields a timescale longer than that of the bulk viscosity 
by a factor of ~ ac s /(/3 Av ). The bulk viscosity therefore dominates the shear for the 
classical AV, provided only that Av << c s . 

If the AV is instead given by HK form [eq. Ql3|)1, dimensional analysis gives 

^«-A^-A^. (HK AV) (35) 



N 



if (V ■ v)j or (V ■ v)j is negative (otherwise 11^ = 0). Although our idealized velocity field 
satisfies (V • v)j = 0, the numerical estimation of the velocity divergence, as computed 
by equation (|T5D , gives small but non-zero results. In deriving equation ( |3l)| ) we have 

1 /2 

used |(vj — Vj) ■ VjWy|/n ~ Av/ (hN N ), which implies |V ■ v|j ~ Av/[hN^ ) from 
eq. flTSp. Before we can estimate vf v = \ — J2j mjHijViWij\ we must note that the 
summation —J2j m j^-ij^iWij appearing in equation fllPf) scales like the number of terms 

1 /2 

Njy in the summation (not as with the classical AV): the condition (V ■ v)j < in 
equation flliD requires that the vectors VjWjj for which TI^ ^ are found preferentially 
in the direction of particle z's velocity deviation from the local fluid flow. Therefore, 
v AV m k[ac s Av/(hN^ 2 ) + k' 2 (3Av 2 /(hN N ), and the timescale satisfies 

where j[, j 2 , fc^ and fc 2 are coefficients of order unity. 



Comparing equations (p3| ) and ( |36| ) we see that the timescale due to the bulk viscosity 
is of the same order of magnitude for the classical and HK artificial viscosities; however, the 
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timescale associated with the von Neumann-Richtmyer term is longer in the HK AV by a 

1 /2 

factor NpJ . Since typical 3D calculations have ~ 50 100, the increase in the viscous 

dissipation timescale is substantial whenever von Neumann-Richtmyer viscosity terms are 
significant. 

If the AV is given by Balsara's form [eq. ([16])], we need to estimate the size of fi 
[eq. (0)] before we can estimate 11^. For our assumed velocity field |V x v| = (2 — X)Q. 
Therefore, provided that A is far enough from 2 that the curl of the velocity dominates 
over the other terms in the denominator on the right hand side of equation ( fL8|) , an SPH 
evaluation of fi gives 

|V ' V| ' ' (37) 



fi 



IV x v| 



N l J\2 - A) 



Recalling that (pi/ pj +Pj/p]) ~ ^ c s/(lP)-> we estimate from equation (|16|) that 



,aAvc 




(Balsara AV) 



(38) 



where j" and are coefficients of order unity. Therefore, i 
2k'^Av 2 f 2 /( 1 hN 1 J 2 ), and the viscous timescale is given by 



AV 



2k'[ac s Avf/^hNU 2 ) + 



\AV 



aAvc* 



hN 



1/2 

iV 



hN 



1/2 
N 



>7 



a\ 2 c* 



--4 k"> 



p\ 3 Av 



-\ -i 



wN N (2-X) 1 zuN 3 J 2 (2- A) 2 7. 



(Balsara AV) (39) 



where k'[ and k'^ are also coefficients of order unity. 



To test these simple analytic estimates we computed Ti = v.-J\ — J2j m j^-ij^iWij\ for 
a rapidly and differentially rotating equilibrium configuration. This configuration was 
constructed in three steps: (1) we created an n = 3, T 1 = 5/3 polytrope (pressure profile 
p = Ap 5 / 3 oc p 4 / 3 , and consequently A oc p 1 / 3 ) of radius R and mass M, (2) assigned 
a velocity t>o = 0.5 (in units where G = M = R = 1) in the azimuthal direction to 
all particles, and (3) relaxed to a rotating equilibrium state by means of an artificial 
"drag" force oc vo<p ~ v i on the particles. The resulting rapidly rotating configuration 
(T/|W| ~ 0.11) is in virial equilibrium with a rotation profile close to Q oc . At small 
w, when the particle smoothing kernels overlap with the rotation axis, the finite resolution 
of the SPH scheme cause deviations from the Q oc w" 1 , cutting off the divergence of fl at 
w = 0. The centrifugal force near w = nevertheless is strong enough to make the density 
a local minimum there; in the equatorial plane the maximum density actually occurs at 
w w 0.14. 
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For such a configuration modeled using N = 10 4 and iVjv ~ 64, Figure 23 compares 
the actual timescale 7$ = — X),- TrijUijV iWij\ computed directly from the SPH code (left 
frame) against our analytic estimates (right frame): (a) classical AV with a = 1, /3 = 0, (b) 
classical AV with a = 0, f3 = 1.5, (c) HK AV with a = 0.5, ,9 = 0, (d) HK AV with a = 0, 
/3 = 0.5, (e) Balsara AV with a = 7/2, p = 0, and (f) Balsara AV with a = 0, = 1.5 x 7/2. 
For all six cases, the same set of particles are analyzed, with the only difference being the 
way vf v is calculated. It is clear that our analytic estimates do a good job of reproducing 
the overall trend in position and magnitude of the actual timescale r. The estimates for 
cases (a) and (c) are identical, while the average measured timescale in case (a) is slightly 
less than that of case (c), which implies k[ < k\. For each of the AV forms, the timescale 
due to the bulk viscosity is significantly less than that due to the von Neumann-Richtmyer 
viscosity. 

Our analytic estimates of and the viscous dissipation timescale r have neglected 
the effects of additional velocity contributions due to particle noise. For this reason, the 
numerical coefficients in equations (0), (0), and (|39"D are not strictly constant but instead 
have some complicated dependence on the neighbor number and noise level in the 
system. Consequently when the particle noise is comparable to the shear velocity, our 
estimates tend to over estimate the timescale. Figure |24] shows the timescales in 6 different 
calculations which have evolved freely for 1 time unit from the relaxed particle state of 
Figure |23|. During this evolution, the particle noise level grows large enough to make our 
analytic formulae overestimate the timescale for cases (d), (e) and (f) by a factor of ~ 2. 
Furthermore, while both the HK and Balsara AVs continue to have significantly longer 
timescales than the classical AV, the timescale for the Balsara AV is now only slightly larger 
than for the HK AV. 

Figure ^ shows the evolution of the angular momentum profile Q in seven different 
calculations which began with the same initial conditions but implemented the different 
artificial viscosities: equations (|TTD, fllHQ and flTBp. The Balsara AV best preserves the 
angular velocity profile. 

One might worry that the spurious increase in the internal energy u or entropy variable 
A due to shear might also occur on as short a timescale as the viscous dissipation. However, 
dimensional analysis on equation ( p0|) and (|2~ID shows that the spurious increase in u and 
A occurs on a timescale ~ TcVild ~~ l)vAv). In typical systems vAv « cf, so that the 
timescale for u or A to change is considerably longer than the viscous dissipation timescale 
t. Figure EB] shows the entropy S as a function of time t for various types of AV. Although 



AVs with more shear viscosity naturally produce more spurious increase in entropy, in all 
cases the rate of entropy increase is rather small. 
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7. Discussion and Summary 
7.1. Particle Diffusion 



Many of our tests focus on spurious diffusion, the motion of SPH particles introduced 
as an artifact of the numerical scheme. Often applications require a careful tracing of the 
particle positions, and in these cases it is essential that spurious diffusion be small. For 
example, SPH calculations can be used to establish the degree of fluid mixing during stellar 
collisions, which is of primary importance in determining the subsequent stellar evolution of 
the merger remnants (e.g. [p5fQ . It must be stressed that the amount of mixing determined 
by SPH calculations is always an upper limit. In particular, low-resolution calculations 
tend to be noisy, and this noise can lead to spurious diffusion of particles, independent of 
any real physical mixing of fluid elements. 

We have analyzed spurious diffusion by using SPH particles in a box with periodic 
boundary conditions to model a stationary fluid of infinite extent. For various noise levels 
(particle velocity dispersions) and neighbor numbers N^, we measure the rate of diffusion, 
which we quantify by calculating a diffusion coefficient D. Although strong shocks and AV 
in SPH calculations can lead to additional particle mixing f28|, particle diffusion is the 



dominant contribution to spurious mixing in weakly shocked fluids. 

Once expressed in terms of the number density of SPH particles and the sound speed, 
these diffusion coefficients can therefore be used to estimate spurious deviations in particle 
positions in a wide variety of applications, including self-gravitating systems. For each 
particle in some large-scale simulation, this spurious deviation is estimated simply from 
equation (|28|) . The coefficient D in the integrand of equation (|28|) depends on the particle's 
velocity deviation from the local flow, the local number density n of particles, and the local 
sound speed c s , so that these quantities need to be monitored for each particle during the 
calculation. Such a scheme is used in §4 to estimate spurious mixing in an equilibrium 
polytrope, and has also been successfully applied in the context of stellar collisions [19| . 



For sufficiently low noise levels, the diffusion coefficient essentially vanishes, as the 
particles simply oscillate around equilibrium lattice sites. We say that such a system has 
"crystallized." For a neighbor number ~ 64, a system of SPH particles will crystallize 
if the root mean square velocity dispersion is less than about 3-4% of the sound speed. We 
find that crystallized cubic lattices are unstable against perturbations, while lattice types 
with large packing fractions, such as hexagonal close-packed, are stable. For this reason 
it may sometimes be better to construct initial data by placing particles in a close-packed 
lattice, rather than in a cubic lattice as is often done. In practice, initial particle data are 
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typically constructed by first relaxing the system with an artificial drag force, a procedure 
which automatically produces a stable lattice structure but also spuriously removes small 
amounts of internal energy. 

The diffusion coefficients have been measured using equal mass particles. Sometimes, 
however, SPH calculations use particles of unequal mass so that less dense regions can still 
be highly resolved. To test the effects of unequal mass particles in a self-gravitating system, 
we constructed an equilibrium n — 1.5 polytrope, using particle masses which increased 
with radius in the initial configuration. Allowing the system to evolve, we observed that 
the heaviest particles gradually migrated towards the center of the star, exchanging places 
with less massive particles. For a polytrope modeled with N « 1.4 x 10 4 particles and a 
neighbor number Njy « 64, the distribution of particle masses is reversed within roughly 
80 dynamical timescales. This is caused by the interactions among neighboring particles 
via the smoothing kernel. These interactions allow energy exchange, and equipartition 
of energy then requires the heavier particles to sink into the gravitational potential well. 
Spurious mixing is therefore a more complicated process in calculations which use unequal 
mass particles: each particle has a preferred direction to migrate, and in a dynamical 
application this direction can be continually changing. For simulations in which fluid mixing 
is important, equal mass particles are an appropriate choice. 



7.2. Shock Tube Tests 



The diffusion tests just described are all done in the absence of shocks and without AV. 
To test the AV schemes described in §2, we turn to a periodic version of the 1-D Riemann 
shock-tube problem of Sod ||56|| . Initially, fluid slabs with constant (and alternating) density 
p and pressure p are separated by an infinite number of planar, parallel, and equally spaced 
interfaces. We treat this inherently 1-D problem with both a 1-D and a 3-D SPH code. 
The 1-D code is naturally more accurate, and provides a benchmark against which we can 
compare the results of our 3-D code. In both cases, periodic boundary conditions allow us 
to model the infinite number of slabs. 

Using various values of a and /3, we performed a number of such shock-tube calculations 
with our 3-D code, at both Mach numbers Ai ~ 1.6 and M. « 13.2 for 7 = 5/3. In 
addition, we performed tests with 7 = 3 and M. ~ 13.2. For each 3-D calculation, we 
compare the time variation of the internal energy and entropy of the system against that 
of the 1-D calculation. Furthermore, since any motion perpendicular to the bulk fluid flow 
is spurious, we were also able to examine spurious mixing. We find that all three forms 
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of AV can handle shocks well. For example, with N = 10 4 and pa 64, there is better 
than 2% agreement with the 1-D code's internal energy vs. time curve when A4 pa 1.6, and 
agreement at about the 3% level when A4 pa 13.2. We also find that both equations (|ITD 
and (|IED , as compared to equation fllBf ), allow less spurious mixing and do somewhat better 
at reproducing the 1-D code's results. For all three forms of AV, increasing the strength of 
the AV allows for less spurious diffusion. 

From Tables 1-4, which present results for numerous shock-tube tests, we see that the 
level at which energy conservation is satisfied depends only weakly on the AV parameters 
but strongly on the length of the timesteps. Energy is typically conserved to better than 
0.1% in the 7 = 5/3 3D calculations whenever the number of timesteps to reach t = 4 
exceeded 1000. Monaghan's timestep routine is more efficient, in part because it takes 
shorter timesteps when shocks are strong (that is, when there are large velocity differentials 
between neighboring particles). The agreement between the 3D and ID calculations for the 
internal energy U and entropy S was strongly dependent on the AV parameters a and f3 
(see §7.4), but only weakly dependent on the Courant number Cat or timestep routine. 

Such calculations are a useful and realistic way to calibrate spurious transport, since 
the test problem, which includes shocks and significant fluid motion, has many of the same 
properties as real astrophysical problems. In fact, the recoil shocks in stellar collisions do 
tend to be nearly planar, so that even the 1-D geometry of the shock fronts is realistic. The 
periodic boundary conditions play the role of gravity in the sense that they prevent the gas 
from expanding to infinity. 



7.3. Shear Flows 

To test the various AV forms in the presence of a shear flow, we impose the so-called 
slipping boundary conditions on a periodic box, as is commonly done in molecular dynamics 
(see, e.g., pBfl ). The resulting "stationary Couette flow" has a velocity field close to 
(v x , v y , v z ) = (voy/L, 0, 0) and allows us to measure the numerical viscosity of the particles. 
As in the shock-tube tests, we also examine spurious mixing in the direction perpendicular 
to the fluid flow. These shear tests therefore allow us to further investigate the accuracy 
of our SPH code as a function of the AV parameters and scheme. We find that both the 
Hernquist & Katz AV [eq. CPL3|) 1 and the Balsara AV [eq. fllfiP l exhibit less viscosity than the 
classical AV [eq. fliTD], While the HK AV produces the smallest numerical viscosity for these 
pure shear flows, it also has the largest spurious diffusion coefficient (see Table IV). The 
product t]D is smallest for the HK AV, indicating that this form is well suited for keeping 
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spurious mixing at a manageable level during calculations involving shear flows. For all 
three forms of the AV, increasing a and /3 tends to damp out the noise and consequently 
decrease spurious mixing, but it also increases the spurious shear viscosity. 

Rotation plays an important role in many hydrodynamic processes. For instance, a 
collision between stars can yield a rapidly and differentially rotating merger remnant. Even 
in the absence of shocks, AV tends to damp away differential rotation due to the relative 
velocity of neighboring particles at slightly different radii, and an initially differentially 
rotating system will tend towards rigid rotation on the viscous dissipation timescale. In 
systems best modeled with a perfect fluid, ideally with a viscous timescale r = oo, any such 
angular momentum transport introduced by the SPH scheme is spurious. 

As a concrete example, we consider an axisymmetric equilibrium configuration 
differentially rotating with an angular velocity profile fl(zu) oc zu~ x , where w is the distance 
from the rotation axis and A is a constant of order unity. We then analytically estimate 
the viscous dissipation timescale for each of the three AVs discussed in §2. These analytic 
estimates are found to closely match numerically measured values of the timescale. Both 
the Hernquist & Katz AV [eq. (|13|)] and the Balsara AV [eq. (|T6D1 yield longer viscous 
timescales than the classical AV [eq. flTT|)1, and hence are better at maintaining the angular 
velocity profile. The Balsara AV does best in this regard. 



7.4. Artificial Viscosity Forms and Parameters 

When choosing values of AV parameters, one must weigh the relative importance of 
shocks, shear, and fluid mixing. For this reason, it is an application-dependent, somewhat 
subjective matter to specify "optimal values" of a and (3. Here, however, we roughly 
delineate the boundaries of the region in parameter space that gives acceptable results. 

Our shock-tube tests of §5 are all done with periodic cells each containing mass 
M = 0.625. We find that the quantity (A(U / E) max ) 2 + ((7 - l)AS max /M) 2 is a convenient 
measure of how well a calculation matches the 1-D code's results for both internal energy 
and entropy (note that (7 — l)AS max /M ~ AA max /A for small AS max ). Values of 
A{U / E) max and AS max are listed in Tables 1 through 4. 

Examination of the final three columns in Table I leads us to the following acceptable 
ranges for a in our 7 = 5/3 low Mach number shock-tube tests: 0.2 ^ a ^ 1 for the 
classical AV, 0.1 <; a ^ 0.5 for the HK AV, and 0.2 ^ 2a/"/ ^ 1 for the Balsara AV. 
If spurious diffusion is not a concern, these ranges for a can all be extended down to a 
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lower limit of a — 0. For a given value of a, the acceptable range of (3 is approximately 
given by 0.8 £ 2a + /3 £ 3.3 for the classical AV, and 0.6 £ 2a + /3 ^ 1.3 for the HK 
AV, and 0.8 <; (2a + (3)2/^ <; 3.3 for the Balsara AV. For parameters in these ranges, all 
three AVs handle the low Mach number shocks with roughly the same level of accuracy. 
When Monaghan's timestep routine is used with Cn = 0.3, values of a and (3 which worked 
particularly well in our low Mach calculations included a = 0.2, (3 — 1 for the classical AV, 
a = 0.3, (3 = 0.5 for the HK AV, and a = 0.5 x 7 /2, f3 = 7/2 for the Balsara AV. 

For our high Mach number tests, inspection of Tables 3 and 4 leads to the following 
acceptable ranges for the AV parameters: 1.3 ^ a + f3 ^ 3.5 for the classical AV, 
1 & & + f3 & 1.6 for the HK AV, and 1.9 & (a + /3)2/ 7 ^ 4 for the Balsara AV. The Balsara 
AV seems capable of handling these high Mach number shocks marginally better than the 
classical AV, and both are more accurate than the HK AV. Values of a and (3 which worked 
particularly well in both of our 7 = 5/3 and 7 = 3 high Mach calculations included a = 1, 
(3 = 1.5 for the classical AV, and a = 2 x 7 /2, (3 = 7 /2 for the Balsara AV. With the HK 
AV, a = 0.5, (3 = 1 worked quite well for 7 = 5/3, as did a = 0.5, (3 = 0.5 for 7 = 3. By 
performing these high Mach calculations for two different values of 7, we have determined 
that the ranges of acceptable AV parameters are only weakly dependent on the equation 
of state for both the classical AV and the HK AV. For the Balsara AV, we find that the 
AV parameters should be scaled with 7, so that softer equations of state require larger AV 
parameters. 

Our shear tests of §6 allow us to further examine the accuracy of our SPH code as a 
function of the AV parameters. Not surprisingly, increasing the strength of the AV tends to 
increase the measured viscosity 77 and decrease the measured spurious diffusion coefficient 
D. The product of the viscosity and the diffusion coefficient provides a convenient (but 
somewhat arbitrary) measure of a calculation's accuracy. We find that increasing a typically 
tends to increase the product r\D in our shear tests, and we consequently choose as the 
"optimal" value of a a relatively small value for which the shock-tube tests (both low and 
high Mach number) give acceptable results. 

The combined results of our shock-tube and shear tests therefore suggest a single set 
of AV parameters which are appropriate in a large number of situations: a ~ 0.5, f3 ~ 1 
for the classical AV, a ~ (3 ~ 0.5 for the Hernquist & Katz AV, and a ~ (3 ~ 7/2 for the 
Balsara AV. We will refer to these parameters as "optimal" ; however, these choices should 
be modified depending on the particular application. For instance, if spurious mixing is 
not a concern and if only weak shocks (M. <^ 2) are expected during a calculation, then a 
smaller value of a is appropriate. Likewise, if strong shocks are expected (M. ^ a few) and 
shear viscosity is not a concern, then a stronger AV is justified. 
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The above recommended values for a and (3 correspond to a somewhat weaker AV 
than is typically suggested in the literature (e.g. a « 1, (3 « 2 for the classical AV). While 
larger AV parameters are appropriate in extreme cases {M. ;> 10), we feel our recommended 
values are slightly more accurate for most situations. Furthermore, since errors do not 
change significantly when the energy equation is integrated instead of the entropy equation 
(the only major difference being a larger AS max for the energy equation, by a roughly 
constant amount, see Table II), we conclude that these "optimal" parameters are insensitive 
to the means by which the thermodynamics is treated. However, we have not tested the 
dependence of the optimal AV parameters on the neighbor number in detail, nor have 
we performed test calculations in which both shear flows and shocks are simultaneously 
occurring. 



Morris & Monaghan |30[] have recently tested the classical AV of equation (|T2J) with a 



time-varying viscosity parameter a, and with (3 = 2a. The evolution of a is determined for 
each particle by a source and decay equation, causing the AV to be significantly active only 
in the immediate vicinity of a shock. The results of their tests are encouraging, and their 
idea of time-varying AV coefficients could be applied to any AV form. 

Our results concerning the various AV forms can be summarized as follows. We find 
that the AVs defined by equations flllD and fll6D do equally well both in their handling of 
shocks and in their controlling of spurious mixing, and do slightly better than equation 
(13|). Furthermore, both equations (|13"D and (16|) do introduce less numerical viscosity than 



equation (p]) . Since use of equation fllBD, Balsara's form of AV, does indeed significantly 
decrease the amount of shear viscosity without sacrificing accuracy in the treatment of 
shocks, we conclude that it is an appropriate choice for a broad range of problems. This is 
consistent with the successful use of Balsara's AV reported by Navarro & Steinmetz |59| in 
their models of rotating galaxies. 

Balsara's viscosity was constructed to be quite similar to the classical AV in form; 
the main difference is that Balsara's form contains a "switch" which suppresses the AV in 
regions of large vorticity. It is a simple matter to generate more sensitive switches than 
the one in equation fll6"|) . For instance, instead of (/j + fj)/2 one could use fcfj [or more 
generally {fifj) k , with k ^ 1). Alternatively, in place of the form function fa defined by 
equation (|T8|), one could use 

(V ■ v)f 

9i = (V-v)? + (Vxv)? + «- (40) 

As expected from scaling analyses such as in §6.2, the viscous dissipation timescale can be 
increased by adopting more sensitive switches such as these. However, such switches also 
tend to allow a faster rate of spurious particle diffusion. We have performed a handful of 
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tests which suggest that such generalizations of Balsara's AV may also handle shocks well, 
although more tests are necessary. 

J.C.L. is supported in part by NSF AST 93-15375 and by a New York Space Grant 
Fellowship. A.S. is supported in part by the Natural Sciences and Engineering Research 
Council of Canada. F.A.R. is supported in part by NSF Grant AST-9618116 and by a Sloan 
Research Fellowship. S.L.S. is supported in part by NSF Grant AST 96-18524 and NASA 
Grant NAG5-7152 at the University of Illinois at Urbana-Champaign. Some computations 
were performed at the Cornell Theory Center. This work was also partially supported by 
the National Computational Science Alliance under Grant AST970022N and utilized the 
NCSA SGI/CRAY POWER CHALLENGEarray and the NCSA SGI/CRAY Origin2000. 



-36- 
REFERENCES 

1. Harlow (1988). 

2. D. H. Porter and P. R. Woodward, Astrophys. J. Supp. 93, 309 (1994). 

3. L. B. Lucy, Astron. J. 82, 1013 (1977). 

4. R. A. Gingold and J. J. Monaghan, Mon. Not. R. Astron. Soc. 181, 375 (1977). 

5. J. J. Monaghan, Annu. Rev. Astron. Astrophys. 30, 543 (1992). 

6. F. A. Rasio and J. C. Lombardi, to appear in a special issue of the Journal of 

Computational and Applied Mathematics (JCAM) on Computational Astrophysics, 
eds. H. Riffert and K. Werner (Elsevier Press, Amsterdam, 1998); |astro-ph / 9805089| 

7. J. J. Monaghan and J. C. Lattanzio, Astophys. J. 375, 177 (1991). 

8. A. Burket, M. R. Bate, and P. Boddenheimer, Mon. Not. R. Astron. Soc. 289, 497 

(1997). 

9. A. F. Nelson, W. Benz, F. C. Adams, and D. Arnett, in press (1998). 

10. A. P. Boss, A. G. W. Cameron, and W. Benz, Icaraus 97, 134 (1992). 

11. M. Herant and W. Benz, Astophys. J. 387, 294 (1992). 

12. D. Garcia-Senz, E. Bravo, and N. Serichol, Astrophys. J. Supp. 115, 119 (1998). 

13. P. Laguna, W. A. Miller, W. H. Zurek, and M. B. Davies, Astophys. J. Lett. 410, L83 

(1993). 

14. N. Katz, D. H. Weinberg, and L. Hernquist, Astrophys. J. Supp. 105, 19 (1996). 

15. P. R. Shapiro, H. Martel, J. V. Villumsen, J. M. Owen, Astrophys. J. Supp. 103, 269 

(1996). 

16. N. Katz, Astophys. J. 391, 502 (1992). 

17. M. Steinmetz, Mon. Not. R. Astron. Soc. 278, 1005 (1996). 

18. J. C. Lombardi, F. A. Rasio, and S. L. Shapiro, Astrophys. J. Lett. 445, L117 (1995). 

19. J. C. Lombardi, F. A. Rasio, and S. L. Shapiro, Astrophys. J. 468, 797 (1996). 



-37- 



20. F. A. Rasio and S. L. Shapiro, Astrophys. J. 401, 226 (1992). 

21. F. A. Rasio and S. L. Shapiro, Astrophys. J. 438, 887 (1995). 

22. M. B. Davies, W. Benz, T. Piran, and F. K. Thielemann, Astrophys. J. 431, 742 (1994). 

23. J. M. Centrella and S. L. W. McMillan, Astrophys. J. 416, 719 (1993). 

24. X. Zhuge, J. M. Centrella, and S. L. W. McMillan, Phys. Rev. D 50, 6247 (1994). 

25. X. Zhuge, J. M. Centrella, and S. L. W. McMillan, Phys. Rev. D 54, 7261 (1996). 

26. R. A. Gingold and J. J. Monaghan, Mon. Not. R. Astron. Soc. 204, 715 (1983). 

27. L. Hernquist and N. Katz, Astrophys. J. Supp. 70, 419 (1989). 

28. J. J. Monaghan, J. Comput. Phys. 82, 1 (1989). 

29. D. Balsara, J. Comput. Phys. 121, 357 (1995). 

30. J. P. Morris and J. J. Monaghan, J. Comput. Phys. 136, 41 (1997). 

31. A. E. Evrard, Mon. Not. R. Astron. Soc. 235, 911 (1988). 

32. L. Hernquist, Astrophys. J. Supp. 64, 715 (1987). 

33. W. Benz, R. L. Bowers, A. G. W. Cameron, and W. H. Press, Astophys. J. 348, 647 

(1990). 

34. W. Benz and J. G. Hills, Astrophys. J. 323, 614 (1987). 

35. A. Sills, J. C. Lombardi, C. D. Bailyn, P. Demarque, F. A. Rasio, and S. L. Shapiro, 

Astrophys. J. 487, 290 (1997). 

36. M. Steinmetz and E. Miiller, Astron. and Astrophys. 268, 391 (1993). 

37. M. B. Davies, M. Ruffert, W. Benz, E. Miiller, Astron. and Astrophys. 272, 430 (1993). 

38. H. Kang et al., Astrophys. J. 430, 83 (1994). 

39. S. C. Smith, J. L. Houser, and J. M. Centrella, Astrophys. J. 458, 236 (1996). 

40. M. Ruffert, M. Rampp, and H.-Th. Janka, Astron. and Astrophys. 321, 991 (1997). 

41. J. J. Monaghan and J. C. Lattanzio, Astron. and Astrophys. 149, 135 (1985). 

42. W. Benz, W. L. Slattery, A. G. W. Cameron, Icarus 66, 515 (1986). 



-38- 



43. J. J. Monaghan, Comp. Phys. Rep. 3, 71 (1985). 

44. R. Klessen, Mon. Not. R. Astron. Soc. 292, 11 (1997). 

45. H. M. P. Couchman, P. A. Thomas, F. R. Pearce, Astrophys. J. 452, 797 (1995). 

46. R. Dave, J. Dubinski, L. Hernquist, New Astronomy 2, 277 (1997). 

47. R. W. Hockney and J. W. Eastwood, Computer Simulations Using Particles (Bristol: 

Adam Hilger, 1988). 

48. F. A. Rasio, PhD Thesis, Cornell University (1991). 

49. N. H. Wells et al, Comp. Phys. 4, 507 (1990). 

50. L. Hernquist, Astrophys. J. 404, 717 (1993). 

51. R. P. Nelson and J. C. B. Papaloizou, Mon. Not. R. Astron. Soc. 265, 905 (1993). 

52. R. P. Nelson and J. C. B. Papaloizou, Mon. Not. R. Astron. Soc. 270, 1 (1994). 

53. A. Serna, J.-M. Alimi, and J.-P. Chieze, Astrophys. J. 461, 884 (1996). 

54. F. A. Rasio and S. L. Shapiro, Astrophys. J. 377, 559 (1991). 

55. M. P. Allen and D.J. Tildesley, Computer Simulations of Liquids (New York: Oxford 

Univ. Press, 1989). 

56. G. A. Sod, J. Comput. Phys. 27, 1 (1978). 

57. R. Courant and K. O. Friedrichs, Supersonic Flow and Shock Waves (New York: 

Springer, 1976). 

58. T. Naitoh and S. Ono, Physics Letters 57A, 448 (1976). 

59. J. F. Navarro and M. Steinmetz, Astrophys. J. 478, 13 (1997). 



This manuscript was prepared with the AAS IATgX macros v3.0. 



- 39 - 



Table I: Low Mach Number Shock-Tube Cases with 7 = 5/3 













steps 


steps 


AE/E 


5l + &l 

y * 






AV 






rli 
<1Z 




TjO 


TjO 


elTj 


r„-2/3i 
[n J 






routine 




a 
P 


routine 


OjV 


+ — 1 
Z — ± 


+ — A 
Z — 4 


+ — A 
Z — 4 


ell Z — 4 


LA(U / a )max 




None 






S 


0.1 


436 


1664 


0.06% 


115.7 


0.111 


0.14 


C 





0.1 


S 


0.1 


413 


1402 


0.04% 


25.25 


0.044 


0.051 


C 





1 


s 


0.8 


38 


140 


4.56% 


3.16 


0.014 


0.025 


C 





2.5 


s 


0.1 


295 


1121 


0.04% 


1.92 


0.014 


0.018 


C 





10 


s 


0.1 


281 


1078 


0.04% 


0.88 


0.030 


0.030 


C 





100 


s 


0.1 


307 


1072 


0.05% 


0.40 


0.064 


0.064 


c 


0.1 


1 


M 


0.3 


163 


572 


0.13% 


2.51 


0.012 


0.019 


c 


0.2 


0.5 


M 


0.3 


145 


523 


0.25% 


2.43 


0.013 


0.019 


c 


0.2 


1 


M 


0.3 


167 


585 


0.11% 


1.81 


0.010 


0.020 


c 


0.2 


1 


M 


0.8 


63 


218 


1.31% 


1.81 


0.013 


0.011 


c 


0.2 


1.25 


M 


0.3 


175 


612 


0.07% 


1.72 


0.011 


0.020 


c 


0.3 


1 


M 


0.3 


170 


604 


0.09% 


1.54 


0.012 


0.020 


c 


0.5 


1 


M 


0.3 


180 


653 


0.08% 


1.10 


0.017 


0.019 


c 


0.5 


1 


M 


0.8 


68 


245 


0.78% 


1.09 


0.016 


0.015 


c 


1 





S 


0.8 


36 


134 


1.41% 


0.78 


0.021 


0.018 


c 


1 


1 


S 


0.8 


39 


164 


1.25% 


0.76 


0.025 


0.020 


c 


1 


1 


M 


0.8 


81 


307 


0.41% 


0.74 


0.025 


0.023 


c 


1 


1.25 


M 


0.3 


221 


832 


0.03% 


0.76 


0.026 


0.025 


c 


1 


2 


S 


0.8 


42 


171 


0.92% 


0.72 


0.028 


0.025 


c 


2 





S 


0.1 


278 


1063 


0.02% 


0.51 


0.035 


0.034 


c 


2 


1 


s 


0.8 


56 


231 


0.72% 


0.52 


0.040 


0.049 


c 


3 





s 


0.1 


275 


1053 


0.01% 


0.41 


0.043 


0.042 


c 


3 


1 


s 


0.8 


79 


329 


2.18% 


0.40 


0.047 


0.066 


c 


10 





s 


0.1 


265 


1035 


0.24% 


0.11 


0.071 


0.068 


HK 





1.25 


M 


0.3 


116 


449 


0.28% 


7.40 


0.016 


0.015 


HK 


0.1 


0.5 


M 


0.3 


111 


440 


0.40% 


8.97 


0.018 


0.016 


HK 


0.1 


0.5 


M 


0.8 


42 


161 


2.79% 


6.95 


0.014 


0.025 


HK 


0.1 


1 


M 


0.8 


45 


171 


2.00% 


4.64 


0.018 


0.014 


HK 


0.1 


2 


M 


0.8 


52 


191 


0.98% 


3.19 


0.025 


0.025 


HK 


0.2 


0.5 


M 


0.3 


117 


463 


0.317c 


4.63 


0.013 


0.012 


HK 


0.2 


0.75 


M 


0.3 


119 


467 


0.24% 


3.95 


0.016 


0.016 


HK 


0.3 


0.5 


M 


0.3 


125 


493 


0.22% 


3.05 


0.016 


0.017 


HK 


0.4 


0.5 


M 


0.3 


135 


534 


0.15% 


2.45 


0.020 


0.022 


UTS 

HK 


0.5 


0.5 


M 


0.3 


145 


572 


0.11% 


1.97 


0.025 


0.027 


HK 


0.5 


1 


M 


0.3 


148 


579 


0.06% 


1.78 


0.029 


0.032 


HK 


0.5 


1 


M 


0.8 


56 


218 


0.43% 


1.80 


0.030 


0.031 


HK 


1 





M 


0.3 


196 


768 


0.01% 


1.21 


0.039 


0.040 


HK 


1 


1 


M 


0.3 


198 


772 


0.02% 


1.13 


0.044 


0.045 


TUT*' 

rlK 


1 

i 


1 
1 


A/T 


U.o 


/ 


001 


U.oz /o 


1.10 


U.U44 


U.U40 


HK 


1 


1 


s 


0.8 


39 


151 


4.70% 


1.36 


0.043 


0.080 


B 





2.5x7/2 


M 


0.3 


192 


687 


0.05% 


5.11 


0.012 


0.011 


B 


0.2x7/2 


0.5x7/2 


M 


0.3 


144 


534 


0.28% 


6.86 


0.021 


0.018 


B 


0.5x7/2 


IX7/2 


M 


0.3 


173 


637 


0.12% 


1.98 


0.010 


0.019 


B 


1x7/2 


0.75X7/2 


M 


0.3 


206 


800 


0.09% 


1.14 


0.018 


0.019 


B 


1x7/2 


1x7/2 


M 


0.3 


211 


811 


0.07% 


1.13 


0.019 


0.020 


B 


IX7/2 


IX7/2 


M 


0.8 


79 


304 


0.54% 


1.08 


0.020 


0.018 


B 


IX7/2 


1.25X7/2 


M 


0.3 


216 


819 


0.05% 


1.12 


0.020 


0.020 


B 


IX7/2 


2x7/2 


M 


0.3 


305 


1195 


0.05% 


0.74 


0.031 


0.031 


B 


2x7/2 





M 


0.3 


233 


855 


0.02% 


1.07 


0.022 


0.023 


B 


2x7/2 


1x7/2 


M 


0.3 


309 


1212 


0.04% 


0.70 


0.032 


0.031 


B 


2x7/2 


1.25x7/2 


M 


0.3 


311 


1213 


0.03% 


0.71 


0.032 


0.032 
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Table II: Low Mach Number Shock-Tube Cases with 7 = 5/3 
(Classical AV, Simple timestep routine, Cm = 0.1) 











steps 


steps 


AE/E 


*2 + *2 












evolution 


to 


to 


at 


[ n -2/3] 






a 




r? 


equation 


t = 1 


t = 4 


t = 4 


at t = 4 


A{U/E) max 


ASmax 





1 


0.01 


entropy 


313 


1172 


0.04% 


3.9 


0.015 


0.016 





1 


0.01 


energy 


335 


1295 


0.04% 


4.1 


0.016 


0.021 


1 





0.01 


entropy 


285 


1082 


0.01% 


0.8 


0.021 


0.023 


1 





0.01 


energy 


309 


1222 


0.01% 


0.8 


0.022 


0.030 


1 


1 


0.002 


entropy 


283 


1076 


0.01% 


0.7 


0.025 


0.025 


1 


1 


0.002 


energy 


306 


1215 


0.01% 


0.8 


0.024 


0.032 


1 


1 


0.01 


entropy 


283 


1076 


0.01% 


0.8 


0.025 


0.025 


1 


1 


0.01 


energy 


306 


1210 


0.01% 


0.8 


0.024 


0.032 


1 


1 


0.05 


entropy 


283 


1079 


0.01% 


0.8 


0.024 


0.025 


1 


1 


0.05 


energy 


307 


1219 


0.02% 


0.8 


0.024 


0.031 


2 


1 


0.01 


entropy 


278 


1061 


0.02% 


0.5 


0.036 


0.035 


2 


1 


0.01 


energy 


304 


1197 


0.00% 


0.5 


0.035 


0.042 


3 


1 


0.01 


entropy 


275 


1053 


0.01% 


0.4 


0.044 


0.044 


3 


1 


0.01 


energy 


304 


1198 


0.00% 


0.4 


0.044 


0.050 
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Table III: High Mach Number Shock-Tube Cases with 7 = 5/3 













steps 


steps 


AE/E 


52+52 

y z 












CLi 




LO 


10 


■it 


[„-2/3l 
[n J 






routine 




a 
P 


routine 




Z — ± 


+ — A 
Z — 4 


+ — A 
Z — 4 


ell Z — 4 


L±{U / £y )max 


A ^ 


None 






M 


0.3 


97 


376 


0.04% 


300.1 


0.207 


0.85 


C 





1 


M 


0.1 


411 


1512 


0.02% 


6.02 


0.028 


0.11 


c 





1 


M 


0.8 


83 


247 


1.49% 


5.46 


0.026 


0.12 


c 





5 


M 


0.8 


124 


371 


1.82% 


1.39 


0.065 


0.14 


c 


0.1 


1 


M 


0.3 


227 


687 


0.07% 


3.84 


0.024 


0.096 


c 


0.2 


0.5 


M 


0.3 


206 


682 


0.05% 


4.57 


0.027 


0.14 


c 


0.2 


1 


M 


0.3 


231 


715 


0.06% 


2.85 


0.024 


0.089 


c 


0.2 


1.25 


M 


0.3 


238 


730 


0.12% 


2.27 


0.027 


0.084 


c 


0.3 


1 


M 


0.3 


233 


746 


0.06% 


2.13 


0.025 


0.085 


c 


0.3 


1.25 


M 


0.3 


243 


763 


0.12% 


1.76 


0.028 


0.081 


c 


0.5 


1 


M 


0.3 


245 


827 


0.05% 


1.38 


0.027 


0.079 


c 


0.5 


1.25 


M 


0.3 


252 


830 


0.10% 


1.26 


0.031 


0.075 


c 


0.5 


2.5 


M 


0.3 


283 


896 


0.16% 


1.06 


0.046 


0.063 


c 


0.7 


1.5 


M 


0.3 


268 


936 


0.08% 


0.94 


0.037 


0.068 


c 


1 





M 


0.8 


97 


386 


0.71% 


1.16 


0.027 


0.127 


c 


1 


1 


M 


0.8 


106 


389 


0.27% 


0.85 


0.033 


0.076 


c 


1 


1.5 


M 


0.3 


292 


1058 


0.05% 


0.82 


0.042 


0.062 


c 


1 


2 


M 


0.3 


299 


1057 


0.08% 


0.79 


0.048 


0.063 


c 


1 


2 


M 


0.8 


112 


397 


0.02% 


0.80 


0.045 


0.069 


c 


2 


2 


M 


0.8 


146 


557 


0.04% 


0.56 


0.059 


0.089 


c 


5 





M 


0.8 


258 


1039 


0.08% 


0.27 


0.077 


0.12 


HK 





1.25 


M 


0.3 


131 


494 


0.28% 


9.71 


0.053 


0.072 


HK 


0.2 


0.5 


M 


0.3 


144 


555 


0.33% 


13.41 


0.043 


0.086 


HK 


0.5 


0.5 


M 


0.3 


186 


727 


0.11% 


4.04 


0.041 


0.080 


HK 


0.5 


1 


M 


0.3 


180 


698 


0.08% 


2.78 


0.060 


0.066 


HK 


1 





M 


0.3 


249 


976 


0.08% 


2.90 


0.029 


0.11 


rliY 


1 

1 


U.zo 


i\/r 


U.o 




y /y 


U.U4/0 


O Oft 


n n/i ft 
U.U4D 


n nco 


HK 


1 


1 


S 


0.8 


44 


163 


3.20% 


1.62 


0.069 


0.093 


HK 


1 


1 


M 


0.3 


238 


941 


0.02% 


1.57 


0.073 


0.083 


HK 


1 


1 


M 


0.8 


89 


350 


0.08% 


1.55 


0.068 


0.077 


B 





2.5X7/2 


M 


0.3 


279 


834 


0.31% 


8.25 


0.030 


0.077 


B 


0.2x7/2 


0.5X7/2 


M 


0.3 


194 


664 


0.10% 


16.99 


0.055 


0.19 


B 


0.5x7/2 


0.75X7/2 


M 


0.3 


243 


854 


0.02% 


5.35 


0.029 


0.13 


B 


0.5x7/2 


IX7/2 


M 


0.3 


254 


857 


0.02% 


4.35 


0.025 


0.11 


B 


IX7/2 


0.75X7/2 


M 


0.3 


293 


1076 


0.02% 


1.88 


0.024 


0.089 


B 


IX7/2 


IX7/2 


M 


0.3 


300 


1106 


0.02% 


1.57 


0.026 


0.074 


B 


IX7/2 


IX7/2 


M 


0.8 


112 


413 


0.33% 


1.62 


0.024 


0.080 


B 


IX7/2 


1.25X7/2 


M 


0.3 


301 


1077 


0.03% 


1.45 


0.028 


0.068 


B 


IX7/2 


1.5 X7/2 


M 


0.3 


306 


1080 


0.05% 


1.40 


0.031 


0.066 


B 


1x7/2 


2x7/2 


M 


0.3 


316 


1100 


0.09% 


1.29 


0.037 


0.064 


B 


2x7/2 





M 


0.3 


403 


1617 


0.03% 


0.91 


0.030 


0.065 


B 


2x7/2 


1x7/2 


M 


0.3 


405 


1577 


0.00% 


0.79 


0.041 


0.058 


B 


2x7/2 


1.25X7/2 


M 


0.3 


406 


1562 


0.01% 


0.81 


0.043 


0.063 



Table IV: High Mach Number Shock- Tube Cases with 7 = 3 









steps 


steps 


AE/E 


81 + 51 






AV 






to 


to 


at 


[ n -2/3] 






routine 


a 




t = 1 


t = 4 


t = 4 


at t = 4 


A(U/E) max 


ASmax 


C 


0.2 


0.5 


238 


849 


1.10% 


1.90 


0.140 


0.073 


C 


0.28 


0.56 


240 


867 


1.01% 


1.27 


0.114 


0.065 


C 


0.3 


1.0 


248 


877 


0.63% 


1.03 


0.081 


0.049 


c 


0.5 


1.0 


261 


938 


0.49% 


0.79 


0.045 


0.039 


c 


0.5 


1.25 


264 


939 


0.37% 


0.81 


0.037 


0.034 


c 


0.7 


1.5 


284 


1013 


0.26% 


0.68 


0.040 


0.024 


c 


0.9 


1.8 


307 


1106 


0.20% 


0.59 


0.047 


0.022 


c 


1.0 


1.5 


313 


1147 


0.24% 


0.58 


0.049 


0.021 


HK 


0.2 


0.5 


184 


708 


1.30% 


3.48 


0.076 


0.042 


HK 


0.28 


0.28 


223 


880 


0.54% 


1.28 


0.061 


0.026 


HK 


0.5 


0.5 


216 


844 


0.62% 


1.43 


0.048 


0.027 


HK 


0.5 


1.0 


214 


836 


0.48% 


1.21 


0.077 


0.023 


HK 


0.7 


0.5 


243 


955 


0.40% 


1.15 


0.076 


0.025 


HK 


0.9 


0.9 


269 


1063 


0.18% 


0.83 


0.115 


0.037 


B 


0.5X7/2 


I.OX7/2 


271 


1014 


0.79% 


1.35 


0.094 


0.060 


B 


0.56X7/2 


0.56X7/2 


286 


1100 


0.85% 


1.41 


0.103 


0.068 


B 


I.OX7/2 


0.75X7/2 


329 


1269 


0.47% 


0.83 


0.039 


0.042 


B 


I.OX7/2 


I.OX7/2 


326 


1240 


0.45% 


0.82 


0.031 


0.038 


B 


I.OX7/2 


1.25X7/2 


324 


1226 


0.39% 


0.81 


0.033 


0.035 


B 


1.8x7/2 


1.8x7/2 


421 


1610 


0.22% 


0.60 


0.066 


0.018 


B 


2.0x7/2 


I.OX7/2 


446 


1722 


0.24% 


0.56 


0.066 


0.018 
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Table V: N = 1000, N N = 64, v /c s = O.I7- 1 / 2 , 7 = 5/3, dt = 0.01 Shear tests 



AV routine 


a 


P 


(1^1)1/2 


r, [Men 2 / 3 ] 


D [c s n-V3] 


r)D [Mc 2 !! 1 / 3 ] 


None 






0.337 


3.0(1) xlO" 4 


2.85(6) 


8.4(4) xlO~ 4 


C 


0.0 


1.00 


0.020 


1.332(l)xl0~ 3 


4.59(5) xl0~ 3 


6.13(7)xl0- 6 


c 


0.0 


2.50 


0.016 


2.763(5) xl0~ 3 


3.89(4) xl0~ 3 


1.07(1) xl0~ 5 


c 


0.0 


10.00 


0.012 


8.71(3) xl0~ 3 


3.73(7) xl0~ 3 


3.25(7)xl0~ 5 


c 


0.3 


0.50 


0.013 


5.60(4) xl0~ 3 


3.91(5)xl0~ 3 


2.19(3)xl0~ 5 


C (v 2 = 0.002) 


0.3 


1.00 


0.013 


6.16(6) xl0~ 3 


3.51(3) xl0~ 3 


2.16(3)xl0~ 5 


c 


0.3 


1.00 


0.013 


6.05(5) xl0~ 3 


3.53(7) xlO- 3 


2.14(5)xl0~ 5 


C {r, 2 = 0.05) 


0.3 


1.00 


0.013 


5.71(4) xl0~ 3 


3.95(4) xl0~ 3 


2.26(3) xlO~ 5 


C 


0.5 


1.00 


0.012 


9.09(10) xl0~ 3 


3.53(7) xl0~ 3 


3.21(7)xl0~ 5 


c 


0.8 


1.25 


0.012 


1.37(3) x 10^2 


3.78(6) xl0~ 3 


5.2(1) xlO~ 5 


c 


1.0 


1.00 


0.012 


1.64(4)xl0~ 2 


3.68(4) xl0~ 3 


6.0(1) xlO~ 5 


c 


1.0 


1.25 


0.012 


1.66(3)xl0- 2 


3.44(9) xlO- 3 


5.7(2) xlO- 5 


c 


2.0 


0.00 


0.010 


3.1(l)xl0~ 2 


3.7(l)xl0~ 3 


1.12(5)xl0~ 4 


c 


3.0 


0.00 


0.009 


4.8(2)xl0~ 2 


3.57(3)xl0~ 3 


1.71(8)xl0~ 4 


HK 


0.0 


1.25 


0.082 


1.72(5) xl0~ 4 


0.17(4) 


3.0(7) xlO- 5 


HK 


0.0 


10.00 


0.038 


5.08(4) xl0~ 4 


2.1(5)xl0~ 2 


1.1(3) xlO" 6 


HK 


0.1 


0.50 


0.066 


4.15(2)xl0~ 4 


0.11(4) 


4.5(17) xlO~ 5 


HK 


0.5 


0.50 


0.024 


1.34(l)xl0~ 3 


5.4(l)xl0~ 3 


7.3(2)xl0- 6 


HK 


0.5 


1.00 


0.025 


1.39(3) xl0~ 3 


5.32(15)xl0~ 3 


7.4(3)xl0- 6 


HK 


1.0 


1.00 


0.022 


2.64(3) xl0~ 3 


5.05(13)xl0~ 3 


1.33(4) xlO~ 5 


B 


0.0 x 7/2 


1.00 x 7/2 


0.026 


2.72(l)xl0~ 4 


1.16(3)xl0~ 2 


3.13(7)xl0~ 6 


B 


0.0 x 7/2 


2.50 x 7/2 


0.023 


4.53(l)xl0~ 4 


7.0(1) xl0~ 3 


3.20(7)xl0~ 6 


B 


0.0 x 7/2 


10.00 x 7/2 


0.020 


1.055(3)xl0~ 3 


5.8(1) xlO~ 3 


6.1(l)xl0~ 6 


B (ri 2 = 0.002) 


0.5 x 7/2 


0.50 x 7/2 


0.013 


2.33(2)xl0~ 3 


4.0(2) xlO" 3 


9.3(5) xlO" 6 


B 


0.5 x 7/2 


1.00 x 7/2 


0.018 


2.25(l)xl0~ 3 


3.42(35) xl0~ 3 


7.7(8) xl0~ 6 


B 


1.0 x 7/2 


0.00 x 7/2 


0.015 


4.22(3) xl0~ 3 


3.15(17) xlO~ 3 


1.33(7) xlO~ 5 


B 


1.0 x 7/2 


1.00 x 7/2 


0.015 


4.33(2)xl0~ 3 


3.93(5)xl0- 3 


1.70(3)xl0- 5 


B 


1.0 x 7/2 


2.00 x 7/2 


0.014 


4.444(7) xl0~ 3 


4.15(8) xl0~ 3 


1.84(3) xlO~ 5 


B 


2.0 x 7/2 


0.00 x 7/2 


0.013 


8.5(1) xl0~ 3 


3.69(5) xlO~ 3 


3.13(5)xl0- 5 


B 


3.0 x 7/2 


0.00 x 7/2 


0.012 


1.584(25) xl0~ 2 


3.82(6)xl0~ 3 


6.0(1) xlO~ 5 
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Table VI: N = 1000, N N = 64, v /c s = 0.57" 1 / 2 , 7 = 5/3, dt = 0.01 Shear tests 



AV routine 


a 


p 


l v y+ v z \l/2 


r? [Mc 3 n 2 /' A \ 


D [c s n L '°\ 


~~ 7~i r a //" „2„ 1 /3l 
r)JJ [Mc^n ' J 


None 






0.128 


4.8(4) 


xlO" 


'S 


0.38(7) 


1.8(4) xlO -5 


C 


0.2 


0.75 


0.029 


9.3(6) 


xlO" 


-3 


2.2(2)xl0~ 2 


2.1(2)xl0~ 4 


C {rp = 0.002) 


0.3 


1.00 


0.026 


1.5(1) 


xlO" 


-2 


2.0(l)xl0~ 2 


3.0(3) xlO~ 4 


C 


0.3 


1.00 


0.026 


1.5(1) 


xlO" 


-2 


2.0(l)xl0~ 2 


2.9(3) xlO~ 4 


C (ri 2 = 0.05) 


0.3 


1.00 


0.026 


1.3(1) 


xlO" 


-2 


2.1(2)xl0~ 2 


2.9(3)xl0- 4 


C 


0.3 


0.50 


0.028 


1.19(10) 


xlO" 


-2 


2.1(l)xl0~ 2 


2.4(3)xl0~ 4 


C 


0.4 


0.50 


0.026 


1.6(2) 


xlO" 


-2 


1.9(2)xl0~ 2 


3.1(4)xl0~ 4 


C 


0.5 


0.50 


0.024 


2.2(2) 


xlO" 


-2 


2.0(l)xl0~ 2 


4.3(5)xl0- 4 


C 


0.5 


1.00 


0.023 


2.5(2) 


xlO" 


-2 


1.8(l)xl0~ 2 


4.5(6)xl0~ 4 


C 


0.8 


1.25 


0.019 


4.5(5) 


xlO" 


-2 


1.7(l)xl0- 2 


7.5(10)xl0- 4 


C 


1.0 


0.25 


0.019 


5.4(7) 


xlO" 


-2 


1.59(7)xl0~ 2 


8.5(12)xlO~ 4 


HK 


0.0 


1.25 


0.079 


2.66(5) 


xlO" 


-4 


0.15(1) 


4.0(3)xl0~ 5 


HK 


0.0 


10.00 


0.063 


1.65(1) 


xlO" 


-3 


7.3(4)xl0- 2 


1.21(7)xlO- 4 


HK 


0.1 


0.50 


0.069 


3.91(5) 


xlO" 


-4 


0.106(4) 


4.1(2)xl0~ 5 


HK 


0.2 


0.50 


0.062 


6.69(7) 


xlO" 


-4 


7.3(3)xl0~ 2 


4.9(2) xlO~ 5 


HK 


0.2 


0.75 


0.061 


7.11(8) 


xlO" 


-4 


6.8(3)xl0~ 2 


4.9(2) xlO~ 5 


HK 


0.3 


0.50 


0.059 


9.7(2) 


xlO" 


4 


5.8(3)xl0~ 2 


5.6(4) xlO~ 5 


HK 


0.4 


0.50 


0.056 


1.28(3) 


xlO" 


-3 


5.4(5)xl0~ 2 


6.9(7)xl0- 5 


HK 


0.5 


0.50 


0.055 


1.66(6) 


xlO" 


-3 


5.5(4)xl0~ 2 


9.2(8)xl0~ 5 


HK 


0.8 


1.25 


0.052 


2.8(1) 


xlO" 


-3 


4.6(7)xl0~ 2 


1.3(2)xl0~ 4 


HK 


1.0 


0.25 


0.051 


3.7(2) 


xlO" 


-3 


4.4(6)xlO~ 2 


1.6(2)xl0- 4 


B 


0.0 x 7/2 


1.00 x 7/2 


0.054 


5.90(4) 


xlO" 


-4 


5.3(3)xl0~ 2 


3.1(2)xl0~ 5 


B 


0.0 x 7/2 


2.50 x 7/2 


0.045 


1.245(9) 


xlO" 


-3 


3.1(3)xl0~ 2 


3.8(4) xlO~ 5 


B 


0.0 x 7/2 


10.00 x 7/2 


0.036 


4.10(2) 


xlO" 


-3 


2.87(6) xl0~ 2 


1.18(3) xlO~ 4 


B (rj 2 = 0.002) 


0.5 x 7/2 


0.50 x 7/2 


0.036 


3.8(2) 


xHT 


-3 


2.3(3)xl0~ 2 


8.7(10) xlO~ 5 


B (ri 2 = 0.05) 


0.5 x 7/2 


0.50 x 7/2 


0.037 


3.6(2) 


xlO" 


-3 


2.2(3)xl0~ 2 


7.8(13) xlO" 5 


B 


0.5 x 7/2 


1.00 x 7/2 


0.036 


4.1(2) 


xlO" 


-3 


2.4(2)xl0~ 2 


9.6(ll)xl0" 5 


B 


0.8 x 7/2 


1.25 x 7/2 


0.032 


7.1(4) 


xlO" 


-3 


2.3(2)xl0~ 2 


1.6(2)xl0~ 4 


B 


1.0 x 7/2 


0.00 x 7/2 


0.031 


8.9(6) 


xlO- 


-3 


2.2(2)xl0~ 2 


2.0(2)xl0~ 4 


B 


1.0 x 7/2 


0.75 x 7/2 


0.030 


9.2(6) 


xlO" 


-3 


2.1(2)xl0~ 2 


1.9(2)xl0~ 4 


B 


1.0 x 7/2 


1.00 x 7/2 


0.030 


9.6(7) 


xlO" 


-3 


2.0(3)xl0~ 2 


2.0(3) xlO~ 4 


B 


1.0 x 7/2 


1.25 x 7/2 


0.030 


9.5(6) 


xlO" 


-3 


2.2(2)xl0~ 2 


2.1(3)xl0~ 4 


B 


1.0 x 7/2 


2.00 x 7/2 


0.030 


9.9(6) 


xlO" 


-3 


2.2(2)xl0~ 2 


2.2(2)xl0~ 4 


B 


2.0 x 7/2 


0.00 x 7/2 


0.024 


2.5(2) 


X ID- 


-2 


2.0(2)xl0- 2 


4.9(7)xl0- 4 
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Table VII: N = 1000,i> /c s = O.ly- 1 / 2 ^ = 5/3, dt = 0.01 Shear tests 
AV routine a N N {^Zi)l/2 v [Mc s n 2 / :i ] D [csn,- 1 / 3 ] V D [Mc^n 1 / 3 ] 



B O.OX7/2 I.OOX7/2 20 0.060 6.63(7)xl0- 4 7.0(3)xl0- 3 4.7(2)xl0- 6 

B O.OX7/2 I.OOX7/2 32 0.037 2.98(2)xl0~ 4 6.7(2)xl0~ 3 2.00(7)xl0~ 6 

B O.OX7/2 I.OOX7/2 64 0.026 2.72(l)xl0~ 4 1.16(3)xl0~ 2 3.13(7)xl0~ 6 

B I.OX7/2 O.OOX7/2 20 0.027 4.85(3)xl0~ 3 5.5(2)xl0~ 3 2.67(10)xl0~ 5 

B I.OX7/2 O.OOX7/2 48 0.017 4.48(2)xl0~ 3 3.85(8)xl0~ 3 1.72(4)xl0~ 5 

B I.OX7/2 O.OOX7/2 64 0.015 4.22(3)xl0~ 3 3.2(2)xl0~ 3 1.33(7)xl0~ 5 

B I.OX7/2 I.OOX7/2 20 0.026 4.92(4)xl0~ 3 5.16(8)xl0~ 3 2.54(5)xl0~ 5 

B I.OX7/2 I.OOX7/2 64 0.015 4.33(2)xl0~ 3 3.93(5)xl0~ 3 1.70(3)xl0~ 5 
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Fig. 1. — The number of particles N(v x ) in velocity bins of width 0.001c s for the equilibrium 
state in a typical simple box test, where c s is the sound speed. The iV = 23 3 particles 
interacted with Nn ~ 64 neighbors and began in a simple cubic lattice configuration 
with noise artificially introduced at t — 0. The solid line shows the best fit Maxwellian, 
corresponding to v rms = 0.404, once the system has reached equilibrium. Deviations from 
this best fit are consistent with statistical fluctuations. 

Fig. 2. — The mean square deviation 5 2 and slope d5 2 /dt as a function of time after an 
equilibrium particle velocity dispersion v rms = 0.069c s has been reached in a typical simple 
box test with Nn = 48 and no AV. At late times, the mean square deviation S 2 increases 
approximately linearly with time, and we define the diffusion coefficient D as the slope of 
this line. Units are discussed in §3.1. 

Fig. 3. — The diffusion coefficient D as a function of the root mean square velocity dispersion 
v rms for various neighbor numbers Nn, as measured by simple box tests in which the SPH 
particles began on a simple cubic lattice. 

Fig. 4. — The diffusion coefficient D near crystallization. Conventions are as in Fig. |3|. At 
t = 0, the SPH particles began on either a simple cubic lattice (data points connected by 
solid lines) or a face centered cubic lattice (data points connected by dashed lines). In this 
regime, D has an obvious dependence on this system's history. 

Fig. 5.— This sequence of cross-sectional slabs, each of thickness Az = 1.02n -1 / 3 , in 
a periodic box of dimension 19n -1 / 3 x 19n -1 / 3 x 19n -1 / 3 demonstrates the instability of 
a simple cubic lattice, (a) At t = the N = 19 3 equal mass SPH particles, each with 
Nn ~ 32 neighbors, are initially motionless with only minuscule deviations (due to numerical 
roundoff errors) from the unstable equilibrium positions of a simple cubic lattice, (b) 
At t = 190n _1 / 3 c^ 1 the particles are in the process of shifting their positions, (c) By 
t = SSOn^ 1 / 3 ^ 1 the particles have settled into a new, stable lattice structure. 

Fig. 6. — The internal energy U, kinetic energy T, total energy E = U + T, and entropy 
S of the iV = 20 3 equal mass particles interacting with ~ 64 neighbors for a calculation 
without AV (solid curve) and a calculation with AV (dashed curve). In contrast to the 
previous simple box tests, the smoothing lengths hi are allowed to vary. The particles began 
on a simple cubic lattice with a Maxwell-Boltzmann velocity distribution. 
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Fig. 7. — The mean square deviation 5 2 and root mean square velocity v rms as a function 
of time for N = 16 3 equal mass SPH particles with ~ 32 in a typical simple box test. 
Here the AV is given by equation QTTj ) with a = 1, (3 = 2 and rf = 0.01. The particles begin 
in a simple cubic lattice with a Maxwell-Boltzmann velocity distribution. The AV drives 
Vrms to zero, so that the mean square deviation 5 2 approaches a constant and the diffusion 
coefficient D = d5 2 /dt becomes zero. 

Fig. 8. — A cross-sectional slab of thickness Az = 0.6n -1 / 3 of the final particle configuration 
for the simple box test presented in Fig. [7|. There are clear dislocations separating the 
different lattice orientations. The initially noisy system has been quenched, or "frozen," into 
a crystal by the AV so quickly that the SPH particles did not have opportunity to settle into 
a single orientation. 

Fig. 9. — The "predicted" (dashed curve) and actual mean square displacement (solid curve) 
for the innermost 6400 particles in an equilibrium n = 1.5 polytrope of mass M and radius R 
modeled with TV = 13949 equal mass particles and N N rs 64. For the top frame equation (0) 
is used to update particle positions, while in the bottom frame equation (H), the XSPH 
method, is implemented. 

Fig. 10. — The fractional spurious change in total energy AE/E, (v/c s ) rms and the mean 
square diffusion distance S 2 as a function of Nn evaluated at a time t = 25(R 3 /GM) 1 ^ 2 
during calculations of an equilibrium n = 1.5 polytrope. Circular data points correspond 
to a polytrope modeled with N = 30000 particles, while square data points correspond to 
those with N = 13949 particles. 

Fig. 11. — Histogram of the average SPH particle mass (m) in five radial bins for the initial 
configuration (dashed curve) and t = 80( J R 3 / GM ) 1/2 configuration (solid curve) during the 
evolution of an equilibrium n — 1.5 polytrope of mass M and radius R. This calculation 
employs N = 13949 particles with Nn ~ 64, Cn = 0.8, the simple timestep routine, and no 
AV. 

Fig. 12. — Density and velocity profiles in a shock-tube test with Mach number M. w 1.6 as 
given by the quasi- analytic solution (solid curve) and our 1-dimensional SPH code (dotted 
curve) at a time t = 0.15. An adiabatic equation of state is used with 7 = 5/3. 
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Fig. 13.— Entropy S in a shock-tube at early times t, as given by the quasi-analytic 
solution (solid line) and our ID SPH code (dotted curve), for the same calculation presented 
in Fig. 0. 



Fig. 14. — The pressure P, entropy variable A, density p and velocity component v x as 
given by our ID code (solid curve) and by one of our 3D calculations (dots) at the relatively 
late time t = 1, for the same shock-tube test shown in Figures [12| and |13[ The bar in the 
lower left corner of the uppermost frame has a total length of 4(/i), where (h) = 0.058 is the 
average smoothing length in the 3D calculation. 



Fig. 15. — Dependence of the results of shock-tube calculations on the AV parameters a 
and (3 for the classical AV with our 3D SPH code: (a) a — 0; (3 — 1 (short), 2.5 (long), 10 
(dot short), (b) (3 = 0; a = 1 (short), 2 (long), 3 (dot short), 10 (dot long), (c) (3 = 1; a = 
(dot), 1 (short dash), 2 (long dash), 3 (dot dash). In all cases if = 0.01. The solid line in 
the bottom two frames corresponds to our benchmark ID calculation. 



Fig. 16. — The average square displacement perpendicular to the bulk fluid flow Sy + 6%, 
the ratio of internal to total energy U/E, and the entropy S for three 7 = 5/3 shock-tube 
calculations with different forms of AV: the classical AV with a = 0.5, (3 = 1 (long dashed 
curve), the HK AV with a = (3 = 0.5 (short dashed curve), and the Balsara AV with 
a = [3 = 7/2 (dotted curve). The solid curve in the bottom two frames results from our ID 
SPH code. 



Fig. 17. — As Fig. [12], but for a higher Mach number (A4 rs 13.2) shock-tube test. The solid 
line is the quasi-analytic solution, while the dotted line is the result of our ID SPH code. 



Fig. 18.— Dependence of the high Mach number shock-tube calculations on the AV 
parameters a and j3 for the classical AV and our 3D SPH code. The different lines represent 
different values, as follows: a — 0, (3 — 1 (dotted curve); a = 0,(3 = 5 (short dashed curve); 
a — 1,(3 — (long dashed curve); a = 5,(3 = (dot-dash). The solid curve is the result of 
the ID calculation presented in Fig. [17. 
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Fig. 19. — As Fig. [TBI but for our high Mach number shock-tube test. 



Fig. 20. — As Fig. [Tpj, but for the high Mach number shock-tube test with 7 = 3. The solid 
line is for the ID calculation, and the others are the results of the 3D calculations with the 
following AV schemes and parameters: dotted, a = 0.5,/? = 1.0, classical AV; short dash, 
a = 0.9, (3 = 1.8, classical AV; long dash, a = 0.28,(3 = 0.56, classical AV; dot-short dash, 
a = 1.0 x 7/2, (3 = 1.0 x 7/2, Balsara AV; dot-long dash, a = 1.8 x 7/2,/? = 1.8 x 7/2, 
Balsara AV; short dash-long dash, a = 0.56 x 7/2, /3 = 0.56 x 7/2, Balsara AV. 



Fig. 21. — Particle velocities in the ^-direction plotted against their y coordinates. Slipping 
boundary conditions at y = ±| are used to maintain the shear flow. The system was relaxed 
without AV for the first 10 time units towards a configuration with v x = Voy/L, v y = v z = 
(the solid line), and then allowed to evolve with AV for another 10 time units to the state 
shown in this figure. Here vq = 0.1c s 7 _1//2 , and L — 1 is the unit of length. We used 
iV = 1000 particles each with Nn = 64 neighbors on average and the classical AV with 
a = 0.5 and [3=1. 



Fig. 22. — The spurious square displacement in the direction perpendicular to the fluid 
flow, energies, and entropy as a function of time in three calculations of a shear flow using 
N = 1000 and N N = 64 with different forms of AV: the classical AV with a = 0.5, (3 = 1 
(long dashed curve), the HK AV with a = (3 = 0.5 (short dashed curve), and the Balsara 
AV with a = /? = 1 x 7/2 (dotted curve). The system was relaxed for the first 10 time 
units towards a situation with (v x ,v y ,v z ) = (0.1c s 'j~ 1 ^ 2 y / L, 0, 0), while from t = 10 to 50 
the system freely evolves with slipping boundary conditions. 



Fig. 23. — The viscous timescale as a function of the distance w from the rotation axis for 
various artificial viscosities in a system which has been relaxed into a rapidly, differentially 
rotating configuration: (a) a = 2 and (3 = in equation ([TT|), (b) a = and (3 = 2 in 



equation ([TT|) , (c) a = 2 and (3 = in equation (|13|), (d) a = and (3 = 2 in equation (|T3|). 
(e) a = 2 x 7/2 and (3 = in equation (|16D, and (f) a = and (3 = 2 x 7/2 in equation ([TB]) . 
Both the actual timescale Tj = Vi/\ — Y,j mjTIyV 'jWij\ computed directly from the SPH code 
(left frame) and the analytic estimate (right frame) are shown. Estimates are computed from 
equation ( |34"D with k\ = = 1 used as an approximation for (a) and (b), from equation ( p6[ ) 
with k[ = k' 2 = 1 for (c) and (d), and from equation ( p9|) with k" = k% = 1 for (e) and (f). 



-50- 



Fig. 24. — The timescale t% = Vi/\ — J2j 77ijIIy VjWij| computed directly from the SPH code 
for various artificial viscosities after 1 time unit of free evolution. The various AV schemes 
and parameters a and (3 are the same as in Fig. E3. 



Fig. 25. — The angular velocity fl as a function of cylindrical radius w at times (a) t — 0, 
(b) t = 1 and (c) t — 10 in seven different calculations which began with the same initial 
conditions but implemented different artificial viscosities, namely, from top to bottom in (c): 
no AV (solid curve), a = and (3 = 2 x 7/2 in equation ( [T6j ) (short dash - long dash), a = 
and /3 = 2 in equation ([13]) (dot - short dash), a = 2x7/2 and (3 = in equation (|16|) (dot - 
long dash), a = and /3 = 2 in equation (TTTp (short dash), a = 2 and /3 = in equation ( O ) 
(long dash), and a = 2 and /5 = in equation (|Tl| ) (dotted). 



Fig. 26. — Entropy S as a function of time for the seven calculations presented in Fig. 
The various line types are as in Fig. 
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